陳蒙蒙，周紹光，劉文靜，程圓娥

（1.河海大學(xué)測(cè)繪科學(xué)與工程學(xué)院，江蘇 南京 211100）

LiDAR點(diǎn)云與遙感影像結(jié)合提取建筑物輪廓

陳蒙蒙1，周紹光1，劉文靜1，程圓娥1

（1.河海大學(xué)測(cè)繪科學(xué)與工程學(xué)院，江蘇 南京 211100）

結(jié)合點(diǎn)云數(shù)據(jù)與遙感圖像數(shù)據(jù)提取建筑物輪廓提出一種新方法。首先對(duì)LiDAR點(diǎn)云高程?hào)鸥窕眯螒B(tài)學(xué)濾波得到地物的高程圖DSM；用漸進(jìn)形態(tài)學(xué)對(duì)高程圖進(jìn)行濾波，得到地形圖；通過(guò)克里金插值對(duì)地面點(diǎn)進(jìn)行內(nèi)插得到DTM；以DSM-DTM=DHM得到粗糙的高度圖；通過(guò)中值濾波得到完善的高度圖，進(jìn)行基于高度信息和彩色信息的區(qū)域生長(zhǎng)，最后進(jìn)行基于方位梯度數(shù)據(jù)的分割，優(yōu)化建筑物邊緣，有效地實(shí)現(xiàn)了建筑物輪廓提取。

克里金插值；DHM；圖割

目前，國(guó)內(nèi)外學(xué)者提出了多種將LiDAR點(diǎn)云與遙感影像結(jié)合提取建筑物輪廓的方法[1-3]。在此基礎(chǔ)上，本文提出另外一種結(jié)合LiDAR點(diǎn)云和遙感影像的建筑物輪廓提取方法。首先進(jìn)行點(diǎn)云濾波，得到地面點(diǎn)；再根據(jù)點(diǎn)云的高程?hào)鸥窕玫礁叱虉D，地面克里金插值獲得DTM；然后利用DSM-DTM=DHM，得到地面高度模型DHM（digital height model）；用區(qū)域生長(zhǎng)和局部方差的方法，獲得建筑物屋頂?shù)拇笾掳邏K，根據(jù)斑塊，獲得建筑物斑塊的邊界緩沖區(qū)；最后根據(jù)影像和DHM信息，構(gòu)建圖割的能量函數(shù)，分割得到建筑物的精確輪廓。

1 算法流程

建筑物輪廓提取算法如圖1所示。

圖1 建筑物輪廓提取算法

1.1 地面點(diǎn)獲取并利用克里金法獲得DTM

1.1.1 點(diǎn)云柵格化和填充缺失數(shù)據(jù)獲取高度圖DSM

對(duì)于空間離散點(diǎn)云，首先將其水平投影，求出所有數(shù)據(jù)點(diǎn)坐標(biāo)的最大值與最小值。再根據(jù)這些最大最小坐標(biāo)值建立一個(gè)與坐標(biāo)軸平行的最小長(zhǎng)方形包圍投影數(shù)據(jù)，將小長(zhǎng)方形劃分為相同大小的一系列小正方形（稱為柵格），再將高程信息賦給柵格。若有幾個(gè)數(shù)據(jù)同時(shí)落入一個(gè)單元格，那么只記錄最小高程，此格網(wǎng)記為gmin。當(dāng)一些單元格里無(wú)高程，對(duì)于這些單元格，解決的策略是將最鄰近的單元格的值賦給此單元格。當(dāng)存在大面積的缺失數(shù)據(jù)時(shí)，此方法會(huì)引起錯(cuò)誤，解決辦法是先找到數(shù)據(jù)缺失的區(qū)域并用與其鄰近的區(qū)域中的最小高程值填充該區(qū)域。具體步驟為：

1）對(duì)數(shù)據(jù)缺失區(qū)域進(jìn)行定位。首先創(chuàng)建二值區(qū)，用“1”代表有值的單元格，用“0”代表無(wú)值的單元格。然后對(duì)二值圖進(jìn)行形態(tài)學(xué)閉運(yùn)算，使用半徑為r的“圓盤”結(jié)構(gòu)元素：

式中，d表示脈沖速度（單位km的脈沖數(shù)），可以由原始數(shù)據(jù)計(jì)算出；c表示單元格的大小，單位是m。閉運(yùn)算以后，由于脈沖密度低導(dǎo)致的縫隙消失，面積較大的數(shù)據(jù)缺失區(qū)域被發(fā)現(xiàn)。

2）用各個(gè)數(shù)據(jù)缺失區(qū)域輪廓線上最小高程值填充相應(yīng)的缺失數(shù)據(jù)，用膨脹后的圖減去原圖可以得到數(shù)據(jù)缺失區(qū)域的輪廓邊界，膨脹的結(jié)構(gòu)元素是半徑為1的“disk”結(jié)構(gòu)元素。

3）對(duì)于上面修改過(guò)的gmin，反復(fù)使用最近鄰的高程值替代無(wú)值的單元格。對(duì)每個(gè)格網(wǎng)進(jìn)行高程賦值后，得到地物的原始高程圖。

1.1.2 對(duì)高程數(shù)據(jù)濾波得到地形數(shù)據(jù)

數(shù)學(xué)形態(tài)學(xué)的基本運(yùn)算是腐蝕和膨脹，其通過(guò)結(jié)構(gòu)元素的鄰域執(zhí)行運(yùn)算。集合X被結(jié)構(gòu)元素B腐蝕表示如下[4]：

腐蝕是為了獲得高程較鄰域最低的點(diǎn)，膨脹是為了獲得高程較鄰域最高的點(diǎn)。開運(yùn)算是基于膨脹和腐蝕的，先腐蝕再膨脹。開運(yùn)算具有一個(gè)很好的特點(diǎn)，能夠去除比指定的窗口尺寸小的突出物體。

選用圓盤型結(jié)構(gòu)元素B，這是由于圓盤型結(jié)構(gòu)元素更加適合去除植被并保持地形不變。對(duì)高程圖進(jìn)行數(shù)學(xué)形態(tài)學(xué)的開運(yùn)算，濾除地表上的樹和高粗差點(diǎn)，h極小轉(zhuǎn)換能夠檢測(cè)和填充所有比周圍低h m以上的凹坑，以去除和填充低粗差區(qū)域。

接著對(duì)高程圖進(jìn)行漸進(jìn)形態(tài)學(xué)濾波，逐漸增大窗口的尺寸，直到找到所有的建筑區(qū)域，選擇最大窗口大于最大的建筑物，濾波去除建筑物，得到地形圖。

1.1.3 克里金插值和高度圖的生成

對(duì)得到的地形圖進(jìn)行克里金插值生成DTM，根據(jù)DSM－DTM=DHM得到粗糙的高度圖。由于中值濾波能有效抑制噪聲[5]，中值濾波的基本原理是把數(shù)字圖像或數(shù)字序列中一點(diǎn)的值用該點(diǎn)的一個(gè)鄰域中各點(diǎn)值的中值代替，讓周圍的像素值接近的真實(shí)值，從而消除孤立的噪聲點(diǎn)，而且可以保持圖像的邊緣特性，不會(huì)使圖像產(chǎn)生顯著的模糊，因此將粗糙的地物高度數(shù)據(jù)進(jìn)行中值濾波，得到完善的地物高度圖。

1.2 基于高度信息和遙感影像提取建筑物

1.2.1 利用高度和光譜特點(diǎn)確定建筑物候選點(diǎn)

對(duì)整個(gè)實(shí)驗(yàn)區(qū)進(jìn)行分析后可知：①房屋高度不低于5 m；②建筑物在高度上與當(dāng)?shù)貥淠竞芙咏虎劢ㄖ镯敳康墓庾V特點(diǎn)與地面很相似。根據(jù)建筑物的高度和光譜分布范圍確定閾值，可分辨出大部分建筑物點(diǎn)。

除上述3個(gè)特點(diǎn)，樹木與房屋頂部還可以通過(guò)鄰域高度最小方差加以區(qū)分：將圖2所示的4個(gè)濾波窗口的定位點(diǎn)順次對(duì)準(zhǔn)高度圖像的每一像素，計(jì)算每一窗口鄰域內(nèi)高度數(shù)據(jù)的方差，求出4個(gè)方差中最小一個(gè)作為當(dāng)前像素處的均勻性指標(biāo)。房屋頂部（尖頂或平頂）的均勻性非常高，而樹木的高度方差則相對(duì)要大得多[7]。對(duì)高度和光譜閾設(shè)置值后，可區(qū)分出大部分建筑物與樹木。

圖2 4個(gè)方向的濾波窗口(黑色表示定位點(diǎn))

1.2.2 區(qū)域生長(zhǎng)獲取建筑物斑塊

建筑物圖斑邊緣不整，存在毛刺、不規(guī)則粘連塊。首先通過(guò)形態(tài)學(xué)腐蝕去除邊緣的不規(guī)則和小錯(cuò)誤，然后利用區(qū)域生長(zhǎng)法重新生長(zhǎng)出建筑物斑塊。建筑物都有足夠大的面積，腐蝕操作去除了邊緣部分后，中央部分可作為生長(zhǎng)階段的種子點(diǎn)。生長(zhǎng)同時(shí)考慮高度色彩的相似性[8]。為了盡可能保證房屋邊緣接近直線，生長(zhǎng)過(guò)程按“圈”執(zhí)行，同時(shí)生長(zhǎng)出多個(gè)新點(diǎn)的步驟才視為有效。區(qū)域生長(zhǎng)實(shí)現(xiàn)的步驟如下：①對(duì)圖像順序掃描，找到第1個(gè)還沒有歸屬的像素，設(shè)該像素為（X0，Y0）；② 以（X0，Y0）為中心，考慮（X0，Y0）的4鄰域像素（X，Y），如果（X0，Y0）滿足生長(zhǎng)準(zhǔn)則，將（X，Y）與（X0，Y0）合并（在同一區(qū)域內(nèi)），同時(shí)將（X，Y）壓入堆棧；③從堆棧中取出一個(gè)像素，把它當(dāng)作（X0，Y0）返回步驟②；④當(dāng)堆棧為空時(shí)，返回步驟①；⑤重復(fù)步驟①、④，直到圖像中的每個(gè)點(diǎn)都有歸屬時(shí)，生長(zhǎng)結(jié)束。

1.3 基于方位梯度數(shù)據(jù)的圖割

為了使得生長(zhǎng)所得建筑物斑塊的邊緣盡可能平直，與影像內(nèi)部邊緣線靠近，對(duì)建筑物斑塊的外邊緣線進(jìn)行優(yōu)化[5,9]。以其中一個(gè)斑塊E為例介紹優(yōu)化過(guò)程、展示優(yōu)化結(jié)果。

1.3.1 計(jì)算影像的梯度模值及方向

依次利用圖3所示的a、b、c、d 4個(gè)模板對(duì)斑塊進(jìn)行卷積，令卷積結(jié)果分別為im1、im2、im3、im4，計(jì)算出im1、im2、im3、im4的絕對(duì)值。每個(gè)像素處有4個(gè)非負(fù)的灰度梯度值，算出每個(gè)像素點(diǎn)的最大灰度梯度值并記錄其對(duì)應(yīng)方向數(shù)。若像素p處im1中數(shù)據(jù)絕對(duì)值最大，則表明該點(diǎn)處有豎直邊緣；類似若im1，im2，im3，im4中數(shù)據(jù)絕對(duì)值最大，則分別表明該點(diǎn)處有“左下右上”邊緣、水平邊緣或“左上右下”邊緣。

圖3 邊緣檢測(cè)模板

1.3.2 對(duì)建筑物外邊緣線進(jìn)行優(yōu)化

E斑塊如圖4a，計(jì)算其附近的梯度模值，得到斑塊附近的梯度方向分布情況。梯度模的取值范圍為[0，100]；梯度方向共有4個(gè)不同的值，分別為水平、豎直及±45°方向。搜索出斑塊輪廓線，膨脹此輪廓線得到優(yōu)化緩沖區(qū)。用緩沖區(qū)的梯度模及方向數(shù)據(jù)構(gòu)建圖，圖割法解算最小割得到房屋E的最終輪廓線，圖4b為這一優(yōu)化結(jié)果，圖4c為提取輪廓線疊加于E斑塊上的效果。

圖4 以E斑塊為例進(jìn)行邊緣線優(yōu)化

對(duì)其他房屋斑塊作類似處理，總體結(jié)果如圖5。圖6為提取輪廓線疊加于原始影像上的效果。大部分的輪廓邊緣線比較平直，與建筑物的外形基本相似。

圖5 邊緣線優(yōu)化

圖6 邊緣疊加到原圖

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

由一些高層住宅性建筑物組成實(shí)驗(yàn)區(qū)域，實(shí)驗(yàn)區(qū)域基本情況有：①建筑物結(jié)構(gòu)復(fù)雜，且部分區(qū)域建筑物被高大植被環(huán)繞遮蔽。②汽車種類較多，部分小汽車的高度和面積均與超小型建筑物接近。③ 道路紋理與建筑物很相似。④空間分辨率SR為0.20 m。

將本文建筑物提取結(jié)果與采用目視解譯的方法手動(dòng)分割原圖作為參考的分割結(jié)果，即人工在原圖上分別標(biāo)注出目標(biāo)區(qū)域的輪廓（ground truth），運(yùn)用算法將實(shí)驗(yàn)中的分割結(jié)果分別與 Ground Truth 作交集檢測(cè)錯(cuò)分的像素點(diǎn)數(shù)，則錯(cuò)分率為錯(cuò)分像素占原圖總像素的比例，驗(yàn)證本文方法的有效性。圖7表示用本文的方法提取的建筑物輪廓；圖8表示Ground Truth得到的結(jié)果，表1為二種方法測(cè)試的對(duì)比值。

表1 本文方法與人工解譯精度對(duì)比

經(jīng)實(shí)驗(yàn)對(duì)比分析發(fā)現(xiàn)，在多次的重復(fù)實(shí)驗(yàn)中，人工解譯的結(jié)果不穩(wěn)定，不能精確地表現(xiàn)建筑物的形狀信息，提取精度遠(yuǎn)不如本文算法，驗(yàn)證了算法的有效性和可靠性。

圖7 本文提取結(jié)果

圖8 人工解譯提取結(jié)果

3 結(jié) 語(yǔ)

本研究中提出了提取建筑物的算法，實(shí)驗(yàn)證明算法的有效性?？紤]到建筑物的復(fù)雜性，在今后的工作中還需要對(duì)二個(gè)方面進(jìn)行更深入的研究和拓展：

1）在城市建筑物密集區(qū)，LiDAR很容易區(qū)分建筑物和其他地物，但在地物分布稀疏的區(qū)域，其分辨地物的能力并不比影像強(qiáng)，自動(dòng)劃分不同區(qū)域是一個(gè)有待解決的預(yù)處理問(wèn)題。

2）從遙感數(shù)據(jù)中獲取建筑物的粗糙斑塊是成功提取建筑物最關(guān)鍵的一個(gè)步驟，但對(duì)斑塊進(jìn)行邊緣優(yōu)化是另一難題。期望的結(jié)果是與建筑物的實(shí)際輪廓盡可能接近且盡可能平滑。如何保證邊緣線平滑是下一步研究重點(diǎn)。

[1] VU TT,YAMAZAKI F,MATSUOKA M, Matsusuoka M.Multi-scale Solution for Building Extraction from LiDAR and Image Data[J].International Journal of Applied Earth Observation and Geoinformation,2009(11):281-289

[2] TURKER M,KOC-SAN D. Building Extraction from Highresolution Optical Spaceborne Images Using the Integration of Support Vector Maching(SVM) Classification, Hough TransFormation and Perceptual Grouping[J].International Journal of Applied Earth Observation and Geoinformation, 2015(34):58-69

[3] 劉戰(zhàn)杰.一種新的基于區(qū)域生長(zhǎng)的彩色圖像分割算法[J].山東大學(xué)學(xué)報(bào)（理學(xué)版）,2010,45(6):76-80

[4] 閆麗麗,許長(zhǎng)輝,高井祥,等. 基于數(shù)學(xué)形態(tài)學(xué)細(xì)化算子的改進(jìn)Canny算法研究[J].測(cè)繪科學(xué), 2010, 35(2):82-84

[5] 吳加富,樊景峰.一種改進(jìn)的快速中值濾波算法[J].濟(jì)源職業(yè)技術(shù)學(xué)院學(xué)報(bào),2009, 8(1):21-23

[6] LI Y, WU H. Adaptive Building Edge Detection by Combining LiDAR Data and Aerial Images[J].The International Archives of the Photogrammetry,Remote Sensing and Spatial Information Sciences,2008,8(37):197-202

[7] 李勇,吳華意.基于形態(tài)學(xué)梯度的機(jī)載激光掃描數(shù)據(jù)濾波方法[J].遙感學(xué)報(bào),2008,12(4):633-639

[8] 林卉,趙長(zhǎng)勝,舒寧.基于Canny算子的邊緣檢測(cè)及評(píng)價(jià)[J].黑龍江工程學(xué)院學(xué)報(bào),2003,17(2):3-6

[9] 劉松濤,殷福亮.基于圖割的圖像分割方法及其新進(jìn)展[J].自動(dòng)化學(xué)報(bào),2012, 38(6):911-922

P237

B

1672-4623（2017）02-0030-03

10.3969/j.issn.1672-4623.2017.02.009

2016-04-05。

項(xiàng)目來(lái)源：國(guó)家自然科學(xué)基金資助項(xiàng)目（41271420/D010702）。

陳蒙蒙，碩士研究生，研究方向?yàn)閿?shù)字圖像處理。