李 昂，黃 煌，夏 煜，沈定濤，王結(jié)臣,3,4

(1.南京大學(xué) 地理與海洋科學(xué)學(xué)院，南京 210023； 2.長(zhǎng)江科學(xué)院 空間信息技術(shù)應(yīng)用研究所，武漢 430010；3.江蘇省地理信息技術(shù)重點(diǎn)實(shí)驗(yàn)室,南京 210023； 4.江蘇省地理信息資源開(kāi)發(fā)與利用協(xié)同創(chuàng)新中心,南京 210023)

1 研究背景

LiDAR(Light Detection and Ranging)是一種結(jié)合全球定位技術(shù)、慣性導(dǎo)航技術(shù)和激光測(cè)距技術(shù)的測(cè)量系統(tǒng)，也是一種能夠快速、高精度、實(shí)時(shí)獲取地表三維信息的技術(shù)手段[1]。利用LiDAR點(diǎn)云進(jìn)行建筑提取能夠快速獲取大區(qū)域建筑物空間分布以及建筑物面積、高度、體量等屬性信息，并實(shí)現(xiàn)大區(qū)域建筑三維重建，可應(yīng)用于城市規(guī)劃、災(zāi)害損失評(píng)估、社會(huì)經(jīng)濟(jì)調(diào)查等多種領(lǐng)域[2]。

基于LiDAR點(diǎn)云的建筑物提取流程主要分為兩步：點(diǎn)云濾波和地物分類。首先點(diǎn)云濾波用來(lái)分離地物點(diǎn)云和地面點(diǎn)云，然后從地物點(diǎn)云集中分離出建筑物點(diǎn)云和其他地物點(diǎn)云。點(diǎn)云濾波的方法可被歸納為4類：數(shù)學(xué)形態(tài)學(xué)方法[3-4]、曲面約束方法[5]、內(nèi)插方法[6]、坡面方法[7]等。在建筑提取方法研究中，將點(diǎn)云轉(zhuǎn)換到數(shù)字地表模型(Digital Surface Model，DSM)或者歸一化的DSM(Normal Digital Surface Model，NDSM)后根據(jù)圖像特征提取建筑物是比較主流的方向。例如Brunn和Weudner[8]在NDSM中根據(jù)高程突變閾值、法向量閾值和面積閾值剔除植被留下建筑；Elberink和Maas[9]通過(guò)特征值從NDSM中分割出植被，然后利用聚類方法得到建筑腳點(diǎn)；王燕燕等[10]、許慧敏等[11]基于NDSM采用支持向量機(jī)(Support Vector Machines，SVM)模糊分類和深度學(xué)習(xí)分類方法對(duì)建筑物進(jìn)行分類研究。為了獲取更多地物特征來(lái)提高分類精度，研究人員嘗試融合高分辨率衛(wèi)星影像等多源數(shù)據(jù)進(jìn)行建筑物提取[12]。從發(fā)展趨勢(shì)上看，研究從簡(jiǎn)單地剔除建筑物以外的其他要素，發(fā)展到利用多種特征進(jìn)行分類，再到利用多源數(shù)據(jù)特征進(jìn)行分類。建筑物提取的本質(zhì)是要根據(jù)目標(biāo)對(duì)象的最大差異進(jìn)行分組，并找到適合目標(biāo)對(duì)象分類的特征。然而，目前的方法大多針對(duì)城市區(qū)域進(jìn)行建筑物提取[13]，農(nóng)村區(qū)域的研究則較少，主要是因?yàn)檗r(nóng)村建筑物與植被冠層的高度基本相同，使得建立農(nóng)村地區(qū)的地物提取規(guī)則有一定難度。

針對(duì)現(xiàn)有LiDAR點(diǎn)云建筑物提取研究中存在的問(wèn)題，本文構(gòu)建了一種基于LiDAR點(diǎn)云的農(nóng)村區(qū)域建筑物面域提取方法，該方法主要包括點(diǎn)云濾波、圖像分割及特征分類3個(gè)步驟。濾波過(guò)程通過(guò)加入梯度及方向約束提高效率和精度；圖像分割過(guò)程針對(duì)NDSM圖像的高程信息，對(duì)分水嶺分割算法進(jìn)行改進(jìn)；特征分類過(guò)程結(jié)合農(nóng)村建筑的特點(diǎn)，設(shè)計(jì)了相應(yīng)的分類特征并利用最大似然方法進(jìn)行分類。最終以計(jì)算Kappa系數(shù)、生產(chǎn)者精度和用戶精度綜合評(píng)價(jià)建筑物分類精度。

圖1 研究區(qū)位置和試驗(yàn)區(qū)LiDAR點(diǎn)云Fig.1 Location of study area and LiDAR point clouddata of the experimental area

2 研究數(shù)據(jù)與方法

2.1 研究區(qū)及數(shù)據(jù)

如圖1(a)所示，本文選擇湖南省益陽(yáng)市泗湖山鎮(zhèn)作為研究區(qū)，該區(qū)域?yàn)榈湫偷闹袊?guó)平原農(nóng)村，區(qū)域內(nèi)包括房屋、樹(shù)木、灌木、農(nóng)田等多種地物，建筑物的形狀大多為矩形，建筑物或河流周?chē)N有樹(shù)木。

從研究區(qū)選擇800 m×800 m的試驗(yàn)區(qū)，其LiDAR點(diǎn)云圖像如圖1(b)所示。該試驗(yàn)區(qū)具有建筑密集區(qū)域和建筑稀疏區(qū)域，并存在較多建筑植被粘連情況，能夠較好地反映農(nóng)村建筑形態(tài)。此外，該試驗(yàn)區(qū)還存在一定的地形起伏區(qū)域，因此對(duì)該試驗(yàn)區(qū)進(jìn)行建筑物提取研究比較具有代表性。

使用德國(guó)TopoSys公司的機(jī)載激光測(cè)量?jī)x器HARRIER 86i來(lái)獲取試驗(yàn)區(qū)的LiDAR點(diǎn)云，激光掃描測(cè)量?jī)x為Riegl LMS-Q680i，脈沖重復(fù)頻率80～400 kHz，光束發(fā)散角45°，飛行高度約1 300 m，平均點(diǎn)云密度2.055 pts/m2。

2.2 基于LiDAR點(diǎn)云的建筑提取方法

2.2.1 方法流程

基于LiDAR點(diǎn)云的建筑提取方法整體流程如圖2所示，可歸納為如下步驟：

(1)構(gòu)建一種范圍約束的形態(tài)學(xué)濾波方法對(duì)LiDAR點(diǎn)進(jìn)行濾波分離地面點(diǎn)和地物點(diǎn)。

(2)對(duì)分離出的地面點(diǎn)和原始點(diǎn)云插值生成數(shù)字高程模型(Digital Elevation Model，DEM)和數(shù)字地表模型(DSM)，利用DEM和DSM計(jì)算得到歸一化的數(shù)字地表模型(NDSM)。

(3)由于試驗(yàn)區(qū)內(nèi)的地物對(duì)象存在較多的建筑植被對(duì)象粘連情況，因此使用改進(jìn)的標(biāo)記分水嶺算法對(duì)NDSM圖像進(jìn)行分割，使地物對(duì)象分割成為獨(dú)立的對(duì)象。

(4)對(duì)單個(gè)地物對(duì)象計(jì)算多種特征，根據(jù)特征進(jìn)行最大似然分類得到建筑物地物對(duì)象。

圖2 基于LiDAR點(diǎn)云的建筑提取方法流程Fig.2 Flowchart of the building extraction method based on LiDAR point cloud data

2.2.2 梯度及梯度方向范圍約束的形態(tài)學(xué)濾波

進(jìn)行LiDAR點(diǎn)云數(shù)據(jù)處理之前應(yīng)先選擇合適的點(diǎn)云數(shù)據(jù)組織方式，常用的LiDAR點(diǎn)云數(shù)據(jù)組織方式有：規(guī)則格網(wǎng)、K-D樹(shù)、八叉樹(shù)、不規(guī)則三角網(wǎng)等[14]。本文使用規(guī)則格網(wǎng)將整個(gè)LiDAR點(diǎn)云的平面覆蓋區(qū)域劃分為規(guī)則排列、大小固定的方格網(wǎng)，根據(jù)試驗(yàn)區(qū)的點(diǎn)云密度確定格網(wǎng)間距1 m。此外，由于機(jī)載LiDAR在掃描過(guò)程中會(huì)產(chǎn)生一定的粗差低點(diǎn)噪聲和粗差高點(diǎn)噪聲。因此，在濾波及地物分割之前通過(guò)設(shè)置高程閾值和距離閾值剔除粗差點(diǎn)，根據(jù)點(diǎn)云高程分布的統(tǒng)計(jì)情況設(shè)定試驗(yàn)區(qū)的高程高值閾值和低值閾值；然后計(jì)算LiDAR點(diǎn)云3×3格網(wǎng)單元中任一點(diǎn)到其他點(diǎn)的最近距離，若該距離大于格網(wǎng)間距的2倍，則該點(diǎn)作為孤立點(diǎn)剔除[15]。

傳統(tǒng)的形態(tài)學(xué)濾波是利用固定大小的窗口對(duì)所有格網(wǎng)進(jìn)行數(shù)學(xué)形態(tài)學(xué)開(kāi)運(yùn)算，剔除小于窗口的突出地物。但由于該方法需要對(duì)所有數(shù)據(jù)點(diǎn)進(jìn)行開(kāi)運(yùn)算，效率較低，并且如果窗口過(guò)大，則會(huì)過(guò)濾掉小的地形起伏；如果窗口過(guò)小，則不能剔除大的建筑物。據(jù)此，Zhang和Whitman[3]提出改進(jìn)的形態(tài)學(xué)濾波，利用逐漸增大的窗口進(jìn)行開(kāi)運(yùn)算，對(duì)每次開(kāi)運(yùn)算的結(jié)果比較運(yùn)算前后高程差值，如果大于閾值，則為地物，否則為地形。為了提高運(yùn)算效率，Zhou等[4]提出的方法通過(guò)梯度來(lái)檢測(cè)地物邊界，只對(duì)地物邊界一定范圍進(jìn)行形態(tài)學(xué)濾波。但是該方法對(duì)于建筑或者樹(shù)木周邊的微地形仍然有過(guò)濾效果，例如圖3中建筑左右兩側(cè)的地形起伏可能被當(dāng)作地物點(diǎn)過(guò)濾。

圖3 一維LiDAR點(diǎn)云及濾波示意圖Fig.3 One-dimensional LiDAR point cloud andfiltering diagram

考慮到建筑物和植被等地面物體呈現(xiàn)出高程突變的特點(diǎn)，而地形起伏的高程變化較小。因此，本文在地物邊界進(jìn)行梯度檢測(cè)的基礎(chǔ)上加入梯度方向信息，只對(duì)邊界點(diǎn)梯度方向上一定范圍的點(diǎn)進(jìn)行濾波，具體為：使用LiDAR點(diǎn)云膨脹運(yùn)算后的高程值減去原來(lái)的高程值作為該點(diǎn)梯度，如式(1)所示。

g=(f⊕w)-f。

(1)

式中：w表示濾波窗口；f?w表示對(duì)圖像f利用窗口w進(jìn)行膨脹運(yùn)算；g為梯度。

梯度運(yùn)算窗口大小為3×3，與之對(duì)應(yīng)地，梯度方向采用LiDAR點(diǎn)3×3窗口范圍內(nèi)最大值的方向，得到8個(gè)方向。遍歷所有格網(wǎng)計(jì)算梯度和梯度方向矩陣后，如果點(diǎn)P的梯度大于一定的閾值，記錄點(diǎn)P的位置及梯度方向，對(duì)點(diǎn)P最大梯度方向上W×W矩陣進(jìn)行標(biāo)記，其中W為最大濾波窗口大小。對(duì)被標(biāo)記的所有格網(wǎng)利用逐漸增大的濾波窗口進(jìn)行開(kāi)運(yùn)算，計(jì)算開(kāi)運(yùn)算前后像元高程差，并與設(shè)定的高差閾值Td進(jìn)行比較，高于閾值則為地物點(diǎn)，否則為地面點(diǎn)。當(dāng)窗口大于最大濾波窗口W時(shí)停止濾波。高差閾值Td計(jì)算方法如式(2)所示，即

Td=(Hmax-Hmin)×scale 。

(2)

式中：Hmax和Hmin分別表示濾波窗口內(nèi)的最大和最小高程值；scale為比例因子，設(shè)置為0.1[4]。

2.2.3 基于改進(jìn)分水嶺算法的NDSM圖像分割

在將原始點(diǎn)云去噪后生成DSM時(shí)為了保留邊緣梯度信息，采用最近鄰法插值，將與像元中心點(diǎn)最靠近的LiDAR點(diǎn)高程值作為像元值。將原始點(diǎn)云濾波后得到的地面點(diǎn)采用曲面擬合法插值生成DEM。為了去除地形的影響，采用歸一化數(shù)字地表模型(NDSM)作為地物分割的基礎(chǔ)，將DSM與DEM數(shù)據(jù)結(jié)果相減得到的NDSM圖像。NDSM計(jì)算如式(3)，即

[NDSM]=[DSM]-[DEM] 。

(3)

式中：[NDSM]、[DSM]、[DEM]表示各模型對(duì)應(yīng)的矩陣，此處的相減，就是2個(gè)矩陣在進(jìn)行相減，得到一個(gè)新的矩陣就是[NDSM]。

分水嶺變換把圖像作為地形表面，在每個(gè)局部最低點(diǎn)位置打洞并讓水從洞中涌出，從低到高淹沒(méi)整個(gè)地形，當(dāng)不同匯水盆地的水聚合到一起時(shí)，修壩阻止其聚合，最終這些大壩的邊界對(duì)應(yīng)圖像的分割線。由于局部梯度的不規(guī)則性，該方法容易計(jì)算過(guò)多的局部極小值從而造成過(guò)分割現(xiàn)象。改進(jìn)的標(biāo)記分水嶺分割算法通過(guò)一定的規(guī)則指定與區(qū)域相對(duì)應(yīng)的局部極小值形成二值標(biāo)記圖像，利用該二值標(biāo)記圖像進(jìn)行形態(tài)學(xué)極小值標(biāo)定，并將其作為原始梯度圖像的局部極小值，實(shí)現(xiàn)梯度重建，然后在重建的梯度上進(jìn)行分水嶺變換。由于改進(jìn)的算法只以標(biāo)記的局部極小值點(diǎn)作為出水點(diǎn)，因此有效地控制了分割后對(duì)象的數(shù)量。

本文針對(duì)NDSM圖像設(shè)計(jì)了基于自動(dòng)標(biāo)記的分水嶺分割算法，其工作流程如下：

(1)對(duì)NDSM圖像進(jìn)行形態(tài)學(xué)梯度計(jì)算，得到梯度圖像G(f)，即

G(f)=(f⊕w)-(fΘw) 。

(4)

式中：fΘw表示對(duì)圖像f利用窗口w進(jìn)行腐蝕操作，窗口w的大小設(shè)置為3×3。

(2)根據(jù)研究區(qū)建筑物和植被對(duì)象的特點(diǎn)，采用局部極大值檢測(cè)標(biāo)記植被對(duì)象和建筑物人字屋頂，平面檢測(cè)標(biāo)記平頂建筑物。局部極大值檢測(cè)判斷若5×5窗口內(nèi)像元中心點(diǎn)高程值為窗口最大高程值，則進(jìn)行標(biāo)記。平面檢測(cè)計(jì)算5×5窗口范圍內(nèi)像元高程值的方差，若小于平面檢測(cè)的閾值，則將中心點(diǎn)設(shè)為標(biāo)記點(diǎn)，并保證一個(gè)平面中心點(diǎn)5像元范圍內(nèi)只有一個(gè)中心點(diǎn)標(biāo)記。在計(jì)算得到標(biāo)記后，建筑物的屋頂可能存在標(biāo)記點(diǎn)過(guò)多的現(xiàn)象，此時(shí)將多個(gè)標(biāo)記點(diǎn)相互連接，計(jì)算并判斷連線上像元值的方差，若小于設(shè)定的閾值，則合并連線兩端的標(biāo)記點(diǎn)，并取連線中心點(diǎn)作為合并后標(biāo)記點(diǎn)的位置。

(3)對(duì)梯度圖像G(f)使用步驟(2)中計(jì)算得到的標(biāo)記進(jìn)行標(biāo)定，并疊加初始地物NDSM圖像進(jìn)行分水嶺變換。標(biāo)記點(diǎn)作為出水點(diǎn)，設(shè)定筑壩條件為：梯度達(dá)到閾值(條件a)；地物NDSM圖像高程出現(xiàn)上升并達(dá)到閾值或者到達(dá)地物邊界(條件b)；出現(xiàn)匯水盆地聚合(條件c)。當(dāng)條件a和條件b同時(shí)實(shí)現(xiàn)時(shí)或者出現(xiàn)條件c時(shí)進(jìn)行筑壩。通過(guò)梯度圖像和地物圖像雙層約束，可以實(shí)現(xiàn)更準(zhǔn)確的分割結(jié)果。

2.2.4 特征值提取與最大似然分類

根據(jù)研究區(qū)實(shí)際情況，地物對(duì)象可分為植被和建筑物2大類。為了區(qū)分2大類地物，對(duì)經(jīng)過(guò)分水嶺分割后形成的獨(dú)立地物對(duì)象提取高程紋理及幾何特征，利用貝葉斯公式計(jì)算地物對(duì)象在每個(gè)特征下屬于某一類的概率，然后累加所有特征下屬于某一類地物的概率，選擇概率最大的類別作為該對(duì)象的最大似然分類結(jié)果。以植被和建筑物兩類地物為例，任一待分類對(duì)象X在n個(gè)特征下屬于植被的概率Tv可累加計(jì)算,如式(5)所示，屬于建筑物概率計(jì)算方法同理。

(5)

式中：Pv，Pb分別表示植被和建筑物的先驗(yàn)概率；Pfv和Pfb分別表示在特征f中，對(duì)象X屬于植被和建筑的概率。

根據(jù)建筑物和植被對(duì)象的高程特點(diǎn)，提取的高程紋理包括以下3類：

(1)最大高程差(Maximum Height Difference，MHD)，即地物對(duì)象的高程最大值與高程最小值的差值。建筑物對(duì)象的高程差通常較小，而植被及其它地物對(duì)象的高程差下相對(duì)較大。

(2)高程標(biāo)準(zhǔn)差(Height Standard Deviation，HSD)，即地物對(duì)象的所有像元高程值的標(biāo)準(zhǔn)差。建筑物對(duì)象高程差異小，因此高程標(biāo)準(zhǔn)差較小，而植被及其它地物的高程標(biāo)準(zhǔn)差較大。

(3)高程梯度標(biāo)準(zhǔn)差(Height Gradient Standard Deviation，HGSD)，首先利用式(4)計(jì)算地物對(duì)象的形態(tài)學(xué)梯度，然后去除對(duì)象邊緣像元后計(jì)算梯度標(biāo)準(zhǔn)差。對(duì)于建筑平頂，去除邊緣后的高程梯度近似為0；對(duì)于建筑人字頂，梯度近似相等；因此建筑物的高程梯度標(biāo)準(zhǔn)差均較小，而植被像元分布的不規(guī)律導(dǎo)致其高程梯度標(biāo)準(zhǔn)差較大。

提取建筑物和植被對(duì)象的幾何特征如下3類：

(1)對(duì)象面積(Area)，地物對(duì)象的像元數(shù)與單個(gè)像元面積的乘積。不同地物對(duì)象處在不同的面積范圍之內(nèi)，通過(guò)對(duì)象面積的約束可以縮小分類目標(biāo)的范圍。

(2)似矩度(Rec)，地物對(duì)象面積Sobj除以地物對(duì)象最小外接矩形面積Srec。由于農(nóng)村建筑大多近似矩形，所以通過(guò)似矩度可以較好地區(qū)分建筑物和其他地物對(duì)象。

(3)粗糙度(Roughness)，地物對(duì)象的頂表面積和其投影面積的比值,建筑物的屋頂一般比較平整，粗糙度較低，植被對(duì)象頂部較為復(fù)雜，具有更高的粗糙度。粗糙度的計(jì)算公式為

Rn=Ssurface/Sprj。

(6)

其中，

式中：Rn為粗糙度；Ssurface表示地物對(duì)象頂面表面積，通過(guò)坡度slope(i)及像元的面積s計(jì)算所得，坡度取3×3窗口中八方向的平均值；n為地物對(duì)象的像元數(shù)；Sprj表示對(duì)象的投影面積。

圖4 試驗(yàn)區(qū)的NDSM圖像高程差和分割后的地物對(duì)象Fig.4 NDSM of the experimental area and groundobjects after segmentation

3 結(jié)果與分析

試驗(yàn)區(qū)的LiDAR點(diǎn)云經(jīng)過(guò)去噪和濾波處理后插值生成DSM和DEM，兩者相減得到的NDSM圖像如圖4(a)所示,對(duì)NDSM圖像進(jìn)行分割得到地物對(duì)象如圖4(b)所示。

選擇試驗(yàn)區(qū)地物總數(shù)的15%作為訓(xùn)練數(shù)據(jù)，結(jié)合航拍影像目視解譯進(jìn)行人工分類，并針對(duì)訓(xùn)練集計(jì)算6個(gè)特征值(MHD、HSD、HGSD、Area、Rec、Rn)的概率分布后，通過(guò)最大似然分類法得到植被和建筑物分類結(jié)果如圖5所示(為了更直觀展示，將航拍影像疊加到分類結(jié)果)。

圖5 試驗(yàn)區(qū)的建筑物和植被分類結(jié)果Fig.5 Classification result of building and vegetation inthe experimental area

為了定量評(píng)估分類結(jié)果的精度，以像元為單位計(jì)算生產(chǎn)者精度(Producer Accuracy,PA)、用戶精度(User Accuracy,UA)、Kappa系數(shù)。由試驗(yàn)數(shù)據(jù)統(tǒng)計(jì)所得：試驗(yàn)區(qū)的地物對(duì)象總像元數(shù)為35 409個(gè)；真實(shí)柵格為建筑的像元數(shù)為24 599個(gè)，為植被的像元數(shù)為10 810個(gè)；分類柵格為建筑的像元數(shù)為25 601個(gè)，為植被的像元數(shù)為9 808個(gè)；分類值與真實(shí)值相同的柵格像元數(shù)為32 561個(gè)，計(jì)算可得Kappa系數(shù)為0.805 2。此外，真實(shí)值與分類值均為建筑物的像元數(shù)為23 676個(gè)，真實(shí)值與分類值均為植被的像元數(shù)為8 885個(gè)，真實(shí)值為植被但分類為建筑物的像元數(shù)923個(gè)，真實(shí)值為建筑物但分類為植被的像元數(shù)1 925個(gè)，計(jì)算可得建筑物分類的PA和UA分別為92.48%和96.25%，植被地物的分類PA和UA分別為90.59%和82.19%。

根據(jù)Kappa系數(shù)>0.8可知，試驗(yàn)區(qū)的建筑物分類取得了較為不錯(cuò)的分類精度。比較建筑物和植被分類的PA和UA，建筑物的分類精度相較植被更高，即建筑物的分類結(jié)果更為準(zhǔn)確。在錯(cuò)誤分類的結(jié)果中，植被誤分類為建筑較少，而建筑誤分類為植被比例相對(duì)較高。這可能是由于濾波過(guò)程或者分割過(guò)程中對(duì)于建筑邊緣的處理不夠完善，造成建筑物邊界像元差異較大，從而影響了特征值的計(jì)算。在后續(xù)研究中，將把重點(diǎn)放在對(duì)象邊界的處理上，以求達(dá)到更好的建筑提取效果。

4 結(jié) 論

本文針對(duì)農(nóng)村地區(qū)地形和地物的特點(diǎn)，構(gòu)建了一套完善的LiDAR點(diǎn)云建筑物提取方法，該方法可用于植被和建筑物高度相近且粘連較多的試驗(yàn)區(qū)。具體結(jié)論總結(jié)如下：

(1)由于試驗(yàn)區(qū)地形高低起伏的特點(diǎn)，將梯度及梯度方向與形態(tài)學(xué)濾波結(jié)合可提高濾波算法的效率和精度。

(2)針對(duì)試驗(yàn)區(qū)植被與建筑物高度相近且粘連較多的情況，使用高程與梯度雙約束的分水嶺筑壩條件對(duì)標(biāo)記的分水嶺算法改進(jìn)，可較好地將NDSM圖像進(jìn)行地物分割。

(3)根據(jù)試驗(yàn)區(qū)建筑物和植被特征選擇契合的高程紋理和幾何特征，利用最大似然估計(jì)方法對(duì)建筑物和植被對(duì)象進(jìn)行分類，得到建筑物分類的生產(chǎn)者精度和用戶精度都>90%，Kappa系數(shù)>0.8的分類效果。