杜守基：融合LiDAR點雲與正射影像的建築物圖割優化提取方法

《測繪學報》

構建與學術的橋樑 拉近與權威的距離

融合LiDAR點雲與正射影像的建築物圖割優化提取方法

杜守基1, 鄒崢嶸1, 張雲生1, 何雪1, 王競雪2

1. 中南大學地球科學與信息物理學院, 湖南 長沙 410083;

2. 遼寧工程技術大學測繪與地理科學學院, 遼寧 阜新 123000

收稿日期：2016-10-25；修回日期：2017-07-12

摘要：提出一種基於圖割演算法的建築物LiDAR點雲與正射影像融合提取方法。首先，利用LiDAR點雲計算3個幾何特徵：平整度、法向量分布和高程紋理一致性。同時利用航空正射影像計算顏色特徵——歸一化植被指數（NDVI）。然後將兩類特徵聯合構建能量函數數據項，綜合數字表面模型（DSM）和NDVI構建平滑項，採用圖割演算法優化得到初始的建築物區域。最後利用初始建築物邊緣一定範圍內的正射影像顏色信息，採用前後景分割的思想進一步優化建築物邊緣。應用ISPRS Vaihingen測試數據進行試驗，結果表明本文方法具有較高的建築物提取精度。

關鍵詞：建築物提取激光點雲正射影像圖割演算法建築物邊緣優化

A Building Extraction Method via Graph Cuts Algorithm by Fusion of LiDAR Point Cloud and Orthoimage

DU Shouji1, ZOU Zhengrong1, ZHANG Yunsheng1, HE Xue1, WANG Jingxue2

Abstract: An automatic building extraction method based on graph cuts algorithm fusing LiDAR point cloud and orthoimage is proposed.Firstly, three geometric features are computed from LiDAR points including flatness, distribution of normal vector and GLCM (grey level co-occurrence matrix) homogeneity of normalized height.NDVI is simultaneously calculated from orthoimage.After that, both kinds of features are combined to construct the data term of energy function, then DSM and NDVI is combined to construct smooth term.Thereafter, graph cuts algorithm is applied to obtain the initial building extraction results.Finally, foreground and background segmentation method is employed to optimize the building boundary based on the orthoimage color information in certain range of the initially detected building boundary.ISPRS Vaihingen dataset is used to evaluate the proposed method.The results reveal that the proposed method can obtain high accuracy of the detection building area.

Key words:building extractionLiDAR point cloudorthoimagegraph cuts algorithmbuilding boundary optimization

建築物是城市空間重要的地理信息要素，從遙感數據中自動提取建築物對於城市三維建模、災害評估、數字地圖更新等應用具有重要的作用[1-2]。建築物的自動提取一直是攝影測量與遙感領域一個重要的研究方向，但是由於城市環境複雜，從遙感數據中自動提取建築物仍然面臨許多困難[3]。

近年來，國內外學者對建築物提取進行了廣泛的研究。根據所用數據的不同，可以分為基於影像數據[4-6]、基於LiDAR數據[7-10]，以及聯合影像與LiDAR數據[11-15]3類方法。第1類方法主要利用影像的光譜信息。由於影像廣泛存在的「同物異譜」、「同譜異物」現象，以及陰影、遮擋等現象，單純利用影像提取的建築物有較多的錯誤，難以實用化[16]。利用LiDAR數據提取建築物主要利用點雲的強度、回波以及幾何信息，雖然自動化程度較高，但相較於影像數據，LiDAR數據的空間解析度較低，提取的建築物輪廓不夠精確。

融合LiDAR點雲與影像數據提取建築物，可以充分利用LiDAR數據精確的三維信息和影像的光譜信息以及解析度高的優勢，二者的融合更加有利於自動提取高精度的建築物信息。文獻[11]提出一種面向對象的航空影像與LiDAR數據融合分類方法，但是該方法依賴於高精度的航空影像分割。文獻[12]首先利用光譜信息分離出植被，然後利用高程信息從地面道路和建築物中分離出建築物。該方法單獨利用光譜信息分離植被，利用高程信息分離建築物，沒有將兩種數據進行有效融合。文獻[13]採用支持向量機融合光譜與幾何特徵進行建築物和植被提取，但是監督分類方法需要大量的標記數據。以上方法主要利用了影像的光譜信息來提高建築物提取精度，而忽略了影像的高解析度特徵來優化建築物邊緣信息。文獻[14]利用超高解析度影像在LiDAR輔助下提取建築物輪廓；文獻[15]融合LiDAR與正射影像進行建築物提取和邊緣規則化。兩種方法都是基於影像上的線特徵提取進行邊緣信息規則化，但是線特徵提取通常不完整，易受陰影區域的影響。針對現有方法存在的問題，基於LiDAR點雲與影像數據的特點，本文引入圖割演算法[17]，提出一種融合幾何信息與顏色信息的建築物提取方法，並利用前後景分割的思想來優化建築物邊緣。在本文方法中，一方面通過圖割演算法來融合幾何與顏色信息，考慮鄰域關係以保證提取結果鄰域內的一致性；另一方面，通過前後景分割的方法基於影像的顏色信息來優化建築物邊緣，以利用影像解析度高的優勢。本文方法主要流程如圖 1所示。

1 融合LiDAR點雲與正射影像的建築物提取1.1 數據預處理

1.1.1 生成DSM、DTM、nDSM

儀器誤差或者航空激光點雲獲取過程中存在的飛鳥，使得獲取的點雲中存在粗差點，影響後續的DSM內插以及幾何特徵計算，首先利用基於統計分析的方法剔除粗差點[18]。在此基礎上，採用文獻[19]的方法將點雲分為地面點和非地面點。然後利用全部點雲內插DSM，利用地面點內插DTM，再利用DSM減去DTM獲取nDSM。由於本文採用的圖割優化過程在DSM格網上實現，因此提出一種保邊緣特徵的DSM內插方法，主要包括以下幾步：

(1) 由點雲平均密度估計DSM格網解析度R。

(2) 含有LiDAR點的格網高程採用其包含的LiDAR點的最大高程值。

(3) 缺少LiDAR點的格網高程通過鄰近格網計算。首先將當前格網周圍8鄰域內的有效格網(含有LiDAR點的格網)的高程按照1.5 m間隔統計直方圖。當某一個區間內含有的格網數大於有效格網數的一半時，直接將這一區間格網高程的平均值賦予當前格網，否則利用8鄰域格網中所有的LiDAR點，按照與格網中心距離反距離加權計算當前格網高程。

(4) 如果8鄰域內沒有有效格網，則繼續擴大鄰域範圍，直到內插出有效高程。

為節省後續處理時間，將非地面點雲投影到DSM格網上，生成一張非地面掩膜Mng以指示非地面區域，也就是建築物提取候選區域。在Mng中，「0」指示地面區域，「1」指示非地面區域。

1.1.2 正射影像陰影增強

NDVI常用於區分植被和非植被區域，本文也將NDVI作為區分植被和建築物的一個特徵。但是陰影區域較暗的顏色會對NDVI的計算產生較大影響[20]，對此，首先採用文獻[21]中的方法檢測陰影區域，然後對陰影區域進行增強。具體做法是利用連通分析將陰影區域聚類，然後對每個聚類的對象分別進行陰影增強，即將影像的RGB轉換到HSI顏色空間，然後按照式(1)對S和I通道進行拉伸，再轉換到RGB空間，以實現對陰影區域的增強。

(1)

1.2 特徵計算

本文主要採用了基於點雲的3個幾何特徵：點雲平整度、法向量分布、基於nDSM影像的紋理一致性，以及基於正射影像的顏色特徵NDVI。下面對採用的特徵進行介紹：

1.2.1 點雲平整度

通常情況下，建築物由多個平面組成，而植被區域表面不規則，因此建築物上的LiDAR點雲具有局部平整性。本文採用協方差矩陣分析來計算點雲這一特徵[22]。

令NP=pi|i=1, 2, …,n表示非地面點，pi=(xi,yi,zi)為其中一個採樣點，Np=pj|pj∈NP是pi的k鄰近點集(本文k取經驗值為15)。Np的3×3協方差矩陣C3×3定義如下

(2)

式中，

為Np點集的重心；|Np|是Np點集中點的數量。

通過特徵分解計算協方差矩陣C3×3的特徵值λ、λ1、λ2(0≤λ≤λ1≤λ2)，點pi的平整度fp定義如下

(3)

fp越小，表明pi更可能位於平面上，也就是更可能屬於建築物點。計算完所有非地面點的特徵fp後，每一個DSM格網的平整度取格網內所有LiDAR點的fp均值，對於沒有對應LiDAR點的格網採用1.1節所提的內插方法進行內插。為保證效率，只對Mng中值為「1」的格網進行內插。

1.2.2 點雲法向量分布

點雲平整度fp描述了局部平面特徵，但是當採樣點位於建築物屋脊處或者屋頂面比較複雜時，fp將會更傾向於指示為植被區域。因此本文採用另一個局部幾何特徵——法向量分布。通常情況下，建築物區域局部點雲的法向量會集中於某一個或某幾個方向，而植被區域局部點雲的法向量更發散。本文提出基於法向量方向分布直方圖統計的平方變異係數來描述這一特徵。令Npn是pi的kn鄰近點集，首先計算Npn中所有點的法向量與豎直方向的夾角αV(αV∈[0, 90°])，然後根據αV將所有點統計到不同的直方圖區間。如果pi位於植被區域，則Npn中的點更傾向於均勻分布到不同的直方圖區間，而建築物反之。因此採用此直方圖的平方變異係數描述法向量分布特徵fN

式中

(4)

式中，ndi是直方圖區間數量；ni(i=1, 2, …,ndi)是每一個區間的點的數量；|Npn|是點集Npn中點的數量。fN越小，表明Npn中的點更傾向於均勻分布，即更可能屬於植被區域，反之，則為建築物區域。對於fN，需要設置的參數為鄰域點數量kn以及直方圖區間數量ndi，本文將kn設置為60，ndi設置為6。計算完所有非地面點的特徵fN後，與平整度一樣，內插每一個DSM格網的fN。

1.2.3 基於nDSM影像的紋理特徵

fp和fN是基於點的特徵，通過點的鄰域計算，所以Np和Npn中的點可能來自於不同的地物，比如建築物和與其靠近的植被。此外，當點雲密度較低時，對於面積小的建築物，採樣鄰域點時很有可能會採樣到其他地物的LiDAR點。因此，本文引入基於格網的特徵來補充基於點的特徵。植被區域高程分布複雜，因此nDSM生成的影像(簡稱歸一化高程影像)中表現為豐富的紋理，而建築物區域高程分布單一，高程影像表現紋理單一，因此本文引入灰度共生矩陣一致性來描述這一紋理特徵fTH(見文獻[23])。fTH越小，表明一致性越差，也就是紋理越豐富，指示該格網更可能屬於植被區域。為了保證不同區域的閾值一致，本文將待提取區域的高差歸一化到0~ 60 m，也就是歸一化高程高於60 m的地物直接認為是建築物，這樣對於所有區域，每一個灰度級代表的高程值一致。為了加快計算速度，只對Mng中值為「1」的格網計算fTH。

1.2.4 顏色相關特徵

影像具有豐富的顏色信息，對於建築物與植被具有較好的區分度。本文採用NDVI來區分植被與建築物，定義如下

(5)

式中，IR、R分別為對應顏色波段的值，採用陰影增強後的影像計算。計算完正射影像的NDVI後，將其賦予每個DSM格網。此外，生成一張基於NDVI的植被區域掩膜Mv，將fNDVI>0.3的格網值設為「1」，表示植被區域，其他區域設為「0」表示非植被區域。

以上介紹的4種特徵fp、fN、fTH、fNDVI取值範圍不同，因此利用Logistic方程進行歸一化，定義如下

(6)

式中，x表示歸一化的中心；k控制Logistic方程的深度；k取值對結果影響較小，根據方程曲線的形態將4種特徵分別設置為-35、2.0、0.2、-10；x通過試驗進行確定，即單獨利用某一個特徵進行建築物提取，通過調整x值，提取結果最優的為該特徵的x取值，對於4個特徵，分別設置為0.06、0.8、18、0.12。

1.3 基於圖割演算法的建築物提取

上節介紹的4個特徵僅僅描述一個點或者一個格網的特徵，並沒有考慮上下文關係，因此本文引入圖割演算法來考慮鄰域關係以保證鄰域內的一致性。圖割演算法基本思想是構建一個權重圖，每一條邊的權表示相應的能量值，然後利用最大流/最小割演算法尋找圖的最小割[17]。圖割演算法的目的是通過最小化如式(7)定義的能量函數為每一個結點賦予一個標記lp。

(7)

式中，為數據項；為平滑項；β為平滑項係數；Dp(lp)表示標記lp符合結點p的程度；Vp,q(lp,lq)表示結點p和q的相似程度。

如式(7)所示，能量函數由數據項和平滑項構成，本文建築物提取的結果是將DSM格網分為兩類，一類是建築物(building)，其他為非建築物(non-building)，即lp∈，因此數據項定義如下：

(8)

式中，λ1、λ2分別為幾何和顏色特徵的權重，λ11、λ12、λ13分別為幾何特徵中fp、fN、fTH的權重，且滿足λ1+λ2=1、λ11+λ12+λ13=1。對於Mng中值為「0」的格網，直接賦予標記為「non-building」，值為「1」的格網，根據式(8)設置數據項。

由於DSM格網高程值在局部區域內具有較高的一致性，而在地物邊緣處則有較明顯的差異，尤其是建築物邊緣，因此本文採用格網高程值來構建平滑項。此外，在建築物附近有高程相近的植被的情況下，可能會出現過度傳播現象，易出現植被誤檢為建築物或者建築物被漏檢的問題。因此本文在平滑項中引入基於NDVI指數的懲罰因子δ，平滑項設置如下

(9)

式中，hp和hq分別為相鄰格網的高程值；常數ε保證分母不為0，可由LiDAR數據的高程精度確定，本文設為0.2 m；δ由植被區域掩膜Mv確定是否給予平滑項，當相鄰格網在Mv中的值不同時，給予懲罰，避免分為一類，本文δ取值0.1；平滑項係數β與建築物提取數據(如建築物高度、建築物屋頂複雜度等)有關，本文將其設置為1。

能量函數定義後，採用標準最大流/最小割演算法[24]求解能量函數的最小割，然後確定每一個DSM格網的標記lp。

1.4 初始結果優化

上一節圖割演算法優化後的結果中可能存在一些低矮地物，如灌木、籬笆等，因此進一步採用高程閾值Th和面積閾值Ta去除這些地物，即高程小於Th的格網直接去除，然後採用連通分析將提取結果聚為不同的對象，當對象面積小於Ta時直接去除。本文Th和Ta分別設置為1.5 m和2.5 m2。

2 建築物邊緣優化

經過圖割優化後提取的建築物結果可以視為一張二值圖。在影像信息的輔助下，本文引入前後景分割的思想，利用圖割演算法對建築物的邊緣提取結果進一步優化，能量函數依然採用式(7)，但數據項和平滑項的設置與初始建築物提取階段有所區別。首先，將初始檢測的建築物區域設為前景，非建築物區域設為後景，然後檢測初始建築物區域的邊緣像素，並在邊緣像素附近劃定緩衝區。通過圖割優化，將前景和後景的標記傳播到緩衝區，從而實現建築物邊緣的精確提取。3個區域的數據項分別設置如下

(10)

平滑項由對應正射影像顏色信息構建，定義如下

(11)

式中，R、G、IR分別為對應顏色波段的值；εcol設置為0.1。定義能量函數後，採用標準最大流/最小割演算法[24]解算能量函數對緩衝區域的格網重新標記。從能量函數的數據項可以看出，本文將初始提取的建築物區域以及非建築物區域設置為固定的值，防止在能量最小化過程中標記發生變化，而主要依賴於平滑項將建築物區域和非建築區域的標記傳播到緩衝區，從而實現緩衝區的精確標記。在此基礎上，進一步利用形態學濾波來消除建築物邊緣的鋸齒現象。

3 試驗與分析3.1 試驗數據

為了驗證本文方法，採用如圖 2所示的ISPRS Vaihingen的數據進行測試，數據採集自德國的Vaihingen地區，具有正射影像和激光點雲數據，正射影像如圖 2(a)所示，解析度為0.09 m，點雲平均密度為4 pts/m2，內插的DSM如圖 2(b)所示，覆蓋面積約2.0 km2，擁有超過1800棟建築物。首先採用如圖 2(c)、(d)、(e)所示的3個區域基準數據測試，正射影像黃色線內為測試區域，參考數據由ISPRS官方提供，其中區域1為形狀複雜的密集歷史建築物區域(圖 2(c))，區域2為樹木環繞的高聳居民樓區域(圖 2(d))，區域3為帶有小附屬建築物的居民區(圖 2(e))。整個數據區域沒有官方的參考數據，筆者採用手工標記的形式在正射影像上標記了所有的房屋區域，以用於驗證本文方法的有效性。

3.2 試驗結果

本文中的λ1、λ2以及λ11、λ12、λ13通過試驗確定，分別取λ1=λ2=0.5、λ11=0.25、λ12=0.5、λ13=0.25。對於ISPRS Vaihingen基準數據和全部數據，採用相同的參數設置。

3.2.1 ISPRS Vaihingen基準數據結果

圖 3列舉了區域1邊緣優化前後建築物區域在正射影像上套合的結果。圖中右上角的區域為圖中心黃色方框內區域放大顯示結果。從圖中可以看出，經過優化後，邊緣更準確，消除了鋸齒現象，並在一定程度上恢復了拐角信息。

為了進一步評價本文方法，採用ISPRS官方提供的參考數據進行定量評價。採用的指標包括基於面積(per-area)、目標(per-object)、面積大於50 m2的目標(per-object>50 m2)的完整率(%)，正確率(%)以及質量(%)(見文獻[25])。結果如表 1所示，與ISPRS測試網站上其他方法結果的對比如表 2所示。

表 1基準數據提取結果Tab. 1Building extraction results of benchmark data

表選項

表 2與其他方法的結果對比Tab. 2Comparison results with other methods

表選項

從表 1和表 2可以看出，本文方法具有一定的優勢。對於基於面積的指標，本文方法取得了最優的完整率，質量指標最高。基於目標的指標中，正確率為100%，說明本文方法沒有誤檢目標，且完整率和質量指標也較優異。基於目標面積大於50 m2的指標，全部達到了100%。此外，通過表 1可以看出，本文方法對於3個環境差異較大的區域都取得了較好的提取結果，也說明了本文方法的穩定性。基於面積的評價結果如圖 4所示，黃色為正確提取，紅色為誤檢，藍色為漏檢。從視覺上看，本文方法取得了較好的結果。誤檢主要來自於建築物附近的植被區域，主要是因為這些區域處於陰影之中，即使進行了陰影增強，也難以有效區分。漏檢建築物主要為較低矮的建築物，漏檢主要原因為：有些建築物較矮小，只含有少量激光點而無法被提取，有些建築物部分被植物遮擋而無法提取。此外，一些屋頂結構複雜的建築物以及在地面點濾波時錯分為地面的建築物也被漏檢。本文優化建築物邊緣信息時，採用的是正射影像，因此，正射影像建築物邊緣糾正的效果也對本文方法的結果產生一定影響。

3.2.2 ISPRS Vaihingen全部數據結果

由於ISPRS Vaihingen全區域沒有提供參考數據，因此筆者通過人機交互方式進行了標記。對於Vaihingen全區域測試了基於面積(per-area)、基於目標(per-object)兩個指標，結果如表 3所示，基於面積的評價結果如圖 5所示。

表 3ISPRS Vaihingen全部數據提取結果Tab. 3Building extraction results for whole area of ISPRS Vaihingen

表選項

從表 3和圖 5可以看出，本文方法在全部大範圍數據上也表現出了較好的性能，說明本文方法具有一定的適應性。從圖 5中還可以看出，主要的漏檢建築物面積較小，另外部分建築物屋頂覆蓋了較多的植被也被漏檢。除了部分誤檢是由植被區域導致以外，還有部分誤檢是由於地面點在點雲濾波時被錯分為非地面，這部分區域被誤檢為建築物。該區域右下角為工廠，有大量貨物堆積在空曠地面上，這些貨物也被誤檢為建築物。

本文方法採用C++實現，運行環境為台式機(Intel i5-3470CPU，主頻：3.20 GHz，內存8 GB)，對於3個測試區域的提取時間分別為12、15和14 s，全部區域的提取時間為660 s，其中耗時最多的步驟為點雲幾何特徵計算。總體上看，本文方法具有較好的效率。

4 結論

本文基於圖割演算法，提出融合LiDAR點雲與正射影像的建築物自動提取方法。本文方法充分利用LiDAR點雲與影像提供的幾何以及顏色特徵，通過圖割演算法進行優化，獲得了鄰域平滑一致的結果。綜合DSM和NDVI構建能量函數平滑項，有效防止了建築物與植被區域的標記過度傳播。利用影像解析度高的特點，基於前後景分割的方法進一步優化建築物邊緣，獲得了更加準確的邊緣信息。通過ISPRS Vaihingen地區的數據進行測試，並與其他方法進行比較，證實了本文方法可以有效地提取建築物，精度較高，具有一定的實用推廣價值。進一步顧及建築物邊緣的規則特性提高建築物的檢測精度是下一步的研究內容。

【引文格式】杜守基, 鄒崢嶸, 張雲生, 等. 融合LiDAR點雲與正射影像的建築物圖割優化提取方法[J]. 測繪學報，2018，47(4)：519-527. DOI: 10.11947/j.AGCS.2018.20160534

喜歡這篇文章嗎？立刻分享出去讓更多人知道吧！