谷延超，范東明

(西南交通大學(xué) 地球科學(xué)與環(huán)境工程學(xué)院，四川 成都 611756)

基于地形的LiDAR數(shù)據(jù)缺失區(qū)域填充方法研究

谷延超，范東明

(西南交通大學(xué) 地球科學(xué)與環(huán)境工程學(xué)院，四川 成都 611756)

針對機(jī)載激光點(diǎn)云數(shù)據(jù)柵格化過程中出現(xiàn)的空白區(qū)域進(jìn)行研究，在對鄰近填充和最低點(diǎn)填充分析的基礎(chǔ)上，提出基于地形的填充方法。對缺失數(shù)據(jù)邊界進(jìn)行一維形態(tài)學(xué)濾波得到其邊界各點(diǎn)地形高，利用邊界點(diǎn)構(gòu)建不規(guī)則三角網(wǎng)進(jìn)行地形內(nèi)插填充。實(shí)驗(yàn)結(jié)果表明，地形填充法可適應(yīng)不同原因造成的區(qū)域數(shù)據(jù)缺失，能夠有效保證高程的連續(xù)性并可提高濾波精度。

機(jī)載LiDAR；柵格化；數(shù)據(jù)填充；不規(guī)則離散點(diǎn)；形態(tài)學(xué)濾波

機(jī)載激光雷達(dá)(Airborne Light Detection and Ranging)作為一種主動(dòng)的新型測量技術(shù)，能夠快速獲取大范圍高精度的空間地理信息，得到了廣泛的應(yīng)用[1-3]。LiDAR的原始三維激光點(diǎn)云作為數(shù)字表面模型，不僅包含地面點(diǎn)，還包含樹木、建筑、車輛等地物信息[4]。濾波是從數(shù)字表面模型中提取地面點(diǎn)，進(jìn)而生成數(shù)字地面模型(DEM)，是LiDAR數(shù)據(jù)處理的重要步驟，同樣也是后續(xù)眾多應(yīng)用的前提。

數(shù)字形態(tài)學(xué)作為圖像處理的重要部分，已經(jīng)廣泛應(yīng)用于點(diǎn)云濾波，并取得了較好的效果[1-3 , 5-7]。數(shù)字形態(tài)學(xué)的濾波算法可以是基于原始LiDAR 點(diǎn)云數(shù)據(jù)，也可以是基于柵格數(shù)據(jù)進(jìn)行，實(shí)驗(yàn)表明，后者無論在時(shí)間還是精度上都優(yōu)于前者。由于點(diǎn)云數(shù)據(jù)是不規(guī)則分布的散點(diǎn)，因此在濾波前須對數(shù)據(jù)進(jìn)行柵格化[2,6]。柵格化數(shù)據(jù)是形態(tài)學(xué)濾波算法的基本數(shù)據(jù)，將影響著濾波精度以及數(shù)據(jù)分層渲染的效果等，尤其是大面積數(shù)據(jù)缺失填充方法對濾波精度的影響，但至今未見此方面詳細(xì)的對比分析。本文詳細(xì)介紹了數(shù)據(jù)柵格化的基本流程，在對現(xiàn)有大范圍數(shù)據(jù)缺失的填充方法分析的基礎(chǔ)上提出了一種新的大面積數(shù)據(jù)缺失填充的方法——地形填充法，并對各種填充方法進(jìn)行詳細(xì)的評價(jià)。

1 點(diǎn)云柵格化

LiDAR獲取的數(shù)據(jù)為一系列離散的空間三維坐標(biāo)及反射強(qiáng)度，點(diǎn)云的位置和間隔在空間上都表現(xiàn)為離散性和不規(guī)則性?；谝?guī)則格網(wǎng)的形態(tài)學(xué)濾波可直接利用現(xiàn)有的形態(tài)學(xué)濾波方法進(jìn)行大數(shù)據(jù)的處理[1]，故需要對原始的不規(guī)則離散數(shù)據(jù)進(jìn)行柵格化。

規(guī)則柵格的格網(wǎng)大小應(yīng)根據(jù)點(diǎn)云數(shù)據(jù)的空間間隔確定，一般情況取平均點(diǎn)云間隔以達(dá)到點(diǎn)云數(shù)據(jù)利用最大化。柵格化過程中，如果一個(gè)柵格內(nèi)落入一個(gè)或幾個(gè)激光腳點(diǎn)，把高程最小值記錄在單元格內(nèi)，若柵格內(nèi)無激光腳點(diǎn)，則需對單元格進(jìn)行填充[2,5,6]。在柵格化過程中，需要生成柵格高程影像和柵格掩膜影像。高程影像記錄落入柵格內(nèi)激光腳點(diǎn)的高程最低值；而掩膜影像記錄柵格有無激光腳點(diǎn)落入，存在激光腳點(diǎn)的掩膜圖像柵格值為1，否則為0。

2 空白柵格填充

空白柵格大致可分為兩類：①由數(shù)據(jù)稀疏造成的點(diǎn)狀或者小面積的數(shù)據(jù)缺失(Ⅰ類數(shù)據(jù)缺失)；②由高吸收性地物或者航帶無重疊度造成的大面積數(shù)據(jù)缺失(Ⅱ類數(shù)據(jù)缺失)。

當(dāng)激光腳點(diǎn)的點(diǎn)間距大于影像的空間分辨率時(shí)，會產(chǎn)生Ⅰ類數(shù)據(jù)缺失[2]，可利用限制搜索半徑的鄰近激光點(diǎn)對高程影像和掩膜影像進(jìn)行填充，搜索半徑可根據(jù)設(shè)備的掃描方式以及柵格大小來設(shè)定，通?？蛇x取1～2倍的柵格大小。

產(chǎn)生Ⅱ類數(shù)據(jù)缺失的原因主要可歸為兩類：①相鄰的掃描帶無重疊區(qū)域；②強(qiáng)吸收性的地物對電磁波的吸收(以水體為主)，填充時(shí)都應(yīng)加以考慮。針對Ⅱ類數(shù)據(jù)缺失通常有兩種策略：鄰近填充和最低點(diǎn)填充[5]。

2.1 鄰近填充

鄰近填充利用相鄰的地物性質(zhì)一致進(jìn)行填充，假設(shè)柵格高程影像為I，點(diǎn)云數(shù)據(jù)集為P，若數(shù)據(jù)缺失區(qū)域?yàn)長，則對L內(nèi)的所有單元格逐一進(jìn)行最近點(diǎn)搜索，步驟如下：

1)以空白單元格lij為中心，統(tǒng)計(jì)搜索半徑r內(nèi)的所有激光腳點(diǎn)集N。

2)若N為空，則擴(kuò)大搜索半徑r，重復(fù)1)直至腳點(diǎn)集N為非空為止。

3)計(jì)算lij與鄰近激光腳點(diǎn)集N中所有腳點(diǎn)的距離，距離最小值所對應(yīng)的高程即為待填充點(diǎn)的高程。

2.2 最低點(diǎn)填充

文獻(xiàn)[5]中認(rèn)為：水面往往接近其周邊范圍內(nèi)的最低點(diǎn)，故提出了最低點(diǎn)填充。對消除Ⅰ類數(shù)據(jù)缺失后的掩膜影像進(jìn)行空白區(qū)域探測，統(tǒng)計(jì)所有的空白區(qū)域，依次對Ⅱ類數(shù)據(jù)缺失利用邊界最低點(diǎn)進(jìn)行填充，步驟如下[5]：

1)對掩膜影像進(jìn)行形態(tài)學(xué)腐蝕運(yùn)算得到腐蝕后的掩膜影像，將原掩膜影像與腐蝕后的掩膜影像相減即可得到數(shù)據(jù)缺失的邊界。

2)提取 Ⅱ 類數(shù)據(jù)缺失的邊界中對應(yīng)高程影像I的最小值，賦給高程影像中對應(yīng)的 Ⅱ 類數(shù)據(jù)缺失區(qū)域。

3)依次對影像進(jìn)行1)、2)步驟，直至所有的Ⅱ類數(shù)據(jù)缺失全部填充完成。

圖1顯示邊界提取的過程。圖1(a)為未經(jīng)任何填充的原始掩膜影像；圖1(b)為經(jīng)過Ⅰ類數(shù)據(jù)缺失填充后，殘留的大面積數(shù)據(jù)缺失即Ⅱ類數(shù)據(jù)缺失；圖1(c)為數(shù)據(jù)缺失區(qū)域腐蝕后的掩膜影像；圖1(d)為Ⅱ類數(shù)據(jù)缺失的邊界。

圖1 邊界提取示意圖

通常情況下，針對由航帶無重疊造成的Ⅱ類數(shù)據(jù)缺失進(jìn)行鄰近填補(bǔ)為較優(yōu)，針對水體吸收造成的Ⅱ類數(shù)據(jù)缺失進(jìn)行最低點(diǎn)填補(bǔ)較優(yōu)[5]，然而在現(xiàn)實(shí)數(shù)據(jù)處理過程中，無法通過空白區(qū)域判斷其產(chǎn)生原因，故不能按照上述分類進(jìn)行填充。

2.3 地形填充

鄰近填充會改變地物本身的屬性，例如尺寸等，而在形態(tài)學(xué)濾波中需要設(shè)定的最大濾波窗口與建筑物的尺寸有關(guān)，最大窗口必須大于濾波區(qū)域內(nèi)最大建筑物的尺寸，這樣才能保證在濾波過程中將建筑物過濾[1,5]。同時(shí)鄰近填充可能在數(shù)據(jù)缺失區(qū)域中部形成明顯的高程跳躍，對后續(xù)的濾波處理會造成一定的影響。

水體造成的缺失可分為由于湖泊和河流造成。由湖泊造成的數(shù)據(jù)缺失采用最低點(diǎn)填充方法往往與現(xiàn)實(shí)地形一致；但河流往往具有一定的地形起伏，尤其在大范圍內(nèi)的河流落差較大，此時(shí)利用最低點(diǎn)填充的數(shù)據(jù)并未與現(xiàn)實(shí)地形保持一致。

綜上所述，最低點(diǎn)填充是為了更好地接近水體的高程即數(shù)字高程模型；而鄰近點(diǎn)填充目的是更好地接近周圍地物或地形，注意到如果鄰近點(diǎn)為地形點(diǎn)，則是將填充點(diǎn)更接近數(shù)字高程模型，當(dāng)其為地物時(shí)，也要在后續(xù)的濾波過程中進(jìn)行剔除，以得到相應(yīng)點(diǎn)的地面高程。故無論最低點(diǎn)填充還是鄰近點(diǎn)填充均可理解為將空白區(qū)域填充到相應(yīng)的地形上。

針對 Ⅱ 類數(shù)據(jù)缺失利用區(qū)域邊界，通過一維形態(tài)學(xué)濾波得到其邊界各點(diǎn)地形，然后利用邊界地形基于TIN對空白區(qū)域進(jìn)行填充。地形填充既能統(tǒng)一不同原因引起的 Ⅱ 類數(shù)據(jù)缺失，又能保證填充區(qū)域內(nèi)部和外部地形的連續(xù)性。地形填充的具體流程如下：

1)對掩膜影像進(jìn)行形態(tài)學(xué)腐蝕運(yùn)算得到腐蝕后的掩膜影像，將原掩膜影像與腐蝕后的掩膜影像相減即可得到數(shù)據(jù)缺失的邊界。

2)對提取的邊界進(jìn)行編碼連接。由于提取的邊界點(diǎn)無拓?fù)潢P(guān)系，而一維形態(tài)學(xué)濾波是根據(jù)周邊地形確定各點(diǎn)高程，故將相鄰的邊界點(diǎn)進(jìn)行編碼連接構(gòu)成一維向量，使得相鄰的邊界點(diǎn)在一維向量中同樣緊鄰。

3)邊界點(diǎn)一維形態(tài)學(xué)濾波。根據(jù)實(shí)驗(yàn)區(qū)域的先驗(yàn)知識，確定最大建筑物的尺寸。以最大建筑物的尺寸為濾波結(jié)構(gòu)元的大小，依次得到各邊界點(diǎn)的地形。圖2所示為某一數(shù)據(jù)缺失區(qū)域提取的邊界地形，一維形態(tài)學(xué)濾波較好地提取了邊界各點(diǎn)的地形。

圖2 某數(shù)據(jù)缺失區(qū)域邊界地形提取示意圖

4)邊界點(diǎn)構(gòu)TIN并進(jìn)行內(nèi)插填充。利用邊界點(diǎn)的像素坐標(biāo)構(gòu)建不規(guī)則三角形，對空白區(qū)域的所有柵格進(jìn)行判斷其所屬的三角形，然后利用不規(guī)則三角形頂點(diǎn)的地形高進(jìn)行內(nèi)插求得各空白柵格的地形。

5)依次對影像進(jìn)行1)、2)、3)、4)步驟，直至所有的Ⅱ類數(shù)據(jù)缺失全部填充完成。

3 評價(jià)分析

為了充分對比和分析3種方法的填充效果以及其對濾波結(jié)果的影響，本文采用目視分析和定量分析相結(jié)合，定量分析采用ISPRS提供的精確分類的點(diǎn)云樣本以及其評價(jià)體系[8]。

3.1 目視分析

ISPRS提供的FSite5中包含由水體吸收和航帶無重疊度而造成的 Ⅱ 類數(shù)據(jù)缺失。圖3中格網(wǎng)大小為2 m，圖3(a)顯示 Ⅱ 類數(shù)據(jù)缺失區(qū)域，圖3(b)、圖3(c)、圖3(d)分別為最低點(diǎn)填充、鄰近填充和地形填充的填充影像。區(qū)域1包含航帶無重疊度而造成的 Ⅱ 類數(shù)據(jù)缺失，呈條帶狀長約1600 m，寬約70 m，且條帶兩側(cè)高差較為明顯，最低點(diǎn)填充與邊界地形高程相差很大;鄰近填充在數(shù)據(jù)缺失區(qū)域中線處存在較為明顯的高程跳躍;而地形填充的效果較好，其高程的連續(xù)性以及與周邊地形的吻合性均較為理想。區(qū)域2包含水體吸收而造成的 Ⅱ 類數(shù)據(jù)缺失，長約50 m，寬約20 m，最低點(diǎn)填充和地形填充能較好地對空白區(qū)域進(jìn)行填充，填充部分與周邊地形吻合較好;而鄰近填充則將河流區(qū)域分塊填充，填充視覺效果較差。

圖3 數(shù)據(jù)缺失區(qū)域及其填充影像

綜合上述，利用最低點(diǎn)填充對于水體效果最佳，但對由航帶無重疊造成的數(shù)據(jù)缺失填充效果較差；鄰近填充會在數(shù)據(jù)缺失區(qū)域中央產(chǎn)生高程不連續(xù)，同時(shí)還改變了地物的尺寸；地形填充法可以兼顧兩者的特點(diǎn)，填充效果從整體而言優(yōu)于前兩種方法，其高程的連續(xù)性及與周邊地形的吻合性都達(dá)到了最佳。

3.2 定量分析

在ISPRS提供的15個(gè)樣本中有5個(gè)樣本涉及到大面積數(shù)據(jù)缺失，分別為Sample21、Sample41、Sample51、Sample52、Sample61，其中Sample51、Sample52兩個(gè)樣本中包含由水體吸收和航帶無重疊度而造成的數(shù)據(jù)缺失，其他樣本均只包含航帶無重疊度而造成的數(shù)據(jù)缺失。保證濾波參數(shù)一致，利用形態(tài)學(xué)濾波對3種填充方法獲得的填充影像進(jìn)行濾波，探討不同填充方法對濾波結(jié)果的影響。

濾波質(zhì)量通過三類誤差進(jìn)行定量評價(jià)，其中第1類誤差(Type Ⅰ)是地面點(diǎn)被誤分為地物點(diǎn)的百分比；第2類誤差(TypeⅡ)是地物點(diǎn)被誤分為地面點(diǎn)的百分比；第3類誤差為總誤差，即被錯(cuò)分的點(diǎn)占整個(gè)數(shù)據(jù)的百分比[8-9]。

表1為樣本的三類誤差統(tǒng)計(jì)，最低點(diǎn)填充TypeⅡ較小但TypeⅠ較大，主要是當(dāng)邊界存在較大高差時(shí)，其將邊界周邊內(nèi)幾乎所有的地面點(diǎn)和非地面點(diǎn)全部標(biāo)記非地面點(diǎn)。Sample52 3種填充方法存在較大的差異，圖4中(a)、(b)、(c)分別為其最低點(diǎn)填充、鄰近填充以及地形填充的誤差分布圖。Sample52下半部由水體吸收造成的條帶狀的Ⅱ類數(shù)據(jù)缺失，其上部存在由航帶無重疊造成的Ⅱ類數(shù)據(jù)缺失。從誤差分布圖可得出：3種填充方法在河流沿

岸無明顯區(qū)別，在因航帶無重疊度而造成的Ⅱ類數(shù)據(jù)缺失區(qū)域有明顯的差異，主要是由于最低點(diǎn)填充和鄰近填充產(chǎn)生的高程不連續(xù)造成的。其他樣本地形填充與鄰近填充相比，其TypeⅠ誤差、TypeⅡ誤差和總誤差都保持一致。從所有樣本匯總來看，地形填充效果最優(yōu)，鄰近填充次之，最低點(diǎn)填充較差。綜上分析可得出：柵格數(shù)據(jù)的填充方法會對濾波精度產(chǎn)生影響，若空白區(qū)域兩側(cè)存在較大高差時(shí)，其影響更為明顯。

表1 樣本的三類誤差統(tǒng)計(jì) %

綜上所述，不同填充方法將會對濾波的精度產(chǎn)生影響，填充方法主要在數(shù)據(jù)空白區(qū)域周邊產(chǎn)生影響；3種填充方法對小范圍水體吸收區(qū)域都較為有效；對航帶無重疊造成的數(shù)據(jù)缺失，鄰近填充和地形填充精度相當(dāng)，但在空白區(qū)域兩側(cè)存在較大高差時(shí)地形填充效果優(yōu)于鄰近填充。

圖4 Sample52不同填充方法誤差分布

4 結(jié)束語

本文介紹了點(diǎn)云數(shù)據(jù)柵格化的基本原理，詳細(xì)討論了柵格化過程中存在的空白柵格及其處理方法，在對不同類型數(shù)據(jù)缺失的填充方法分析的基礎(chǔ)上提出了基于地形填充法。利用ISRPS提供的測試數(shù)據(jù)，通過目視分析和定量分析對3種填充方法進(jìn)行了測試和分析，得出：①最低點(diǎn)填充和鄰近填充都會產(chǎn)生高程不連續(xù)現(xiàn)象，而地形填充效果從整體而言優(yōu)于前兩種方法，其高程的連續(xù)性及與周邊地形的吻合性均較理想；②填充方法將影響濾波的精度，填充方法主要在數(shù)據(jù)空白區(qū)域周邊產(chǎn)生影響；③地形填充能夠較好地解決由水體吸收和航帶無重疊等不同原因造成的數(shù)據(jù)缺失，可進(jìn)一步提高濾波精度。

[1]PINGEL T, CLARKE K, MCBRIDE W. An improved simple morphological filter for the terrain classification of airborne LIDAR data [J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2013, 77: 21-30.

[2]沈晶，劉紀(jì)平，林祥國．用形態(tài)學(xué)重建方法進(jìn)行機(jī)載LiDAR 數(shù)據(jù)濾波[J]．武漢大學(xué)學(xué)報(bào):信息科學(xué)版，2011，36 (2)：167-176.

[3]李峰，崔希民，袁德寶，等．改進(jìn)坡度的LiDAR 點(diǎn)云形態(tài)學(xué)濾波算法[J]．大地測量與地球動(dòng)力學(xué)，2012，32 (5)：128-132.

[4]周曉明，馬秋禾，許嘵亮，等．LIDAR點(diǎn)云濾波算法分析—以ISPRS測試實(shí)驗(yàn)為參考 [J]．測繪工程，2011，20 (5)：36-39.

[5]CHEN Q, GONG P, BALDOCCHI D, et al. Filtering airborne laser scanning data with morphological methods[J]. Photogrammetric Engineering and Remote Sensing, 2007, 73(2): 175-184.

[6]羅伊萍，姜挺，王鑫,等．基于數(shù)學(xué)形態(tài)學(xué)的LiDAR 數(shù)據(jù)濾波新方法[J]．測繪通報(bào)，2011(3)：15-19.

[7]隋立春，張熠斌，柳艷，等．基于改進(jìn)的數(shù)學(xué)形態(tài)學(xué)算法的LiDAR點(diǎn)云數(shù)據(jù)濾波[J]．測繪學(xué)報(bào)，2010，39(4)：390-396.

[8]SITHOLE G, VOSSELMAN G. Experimental comparison of filter algorithms for bare-earth extraction from airborne laser scanning point clouds[J]． ISPRS Journal of Photogrammetry and Remote Sensing，2004，59(1-2):85-101.

[9]程曉光，黃先鋒，張帆.機(jī)載LiDAR數(shù)據(jù)的城區(qū)樹木點(diǎn)提取方法[J].測繪科學(xué)，2014，39(3)：52-56.

[責(zé)任編輯：劉文霞]

Research on filling aggregated missing data of LiDAR with topographic method

GU Yan-chao, FAN Dong-ming

(School of Geoscience and Environment Engineering, Southwest Jiaotong University, Chengdu 611756, China)

Through studying on aggregated missing data of LiDAR in the procedures of rasterizing, a method which fills aggregated missing data with terrain of boundary is presented based on analyzing lowest method and nearest method. Terrain of boundary is obtained by one-dimensional morphological filter, and triangulate irregular network (TIN) is constructed. Interpolation is carried out based on Terrain of boundary and TIN. Experiments show that the proposed method adapts to aggregated missing data caused by different reasons to effectively ensure the continuity of the elevation, and improves the filtering accuracy.

airborne LiDAR; rasterizing; data filling; irregular distribution points; morphological filter

2013-08-18,2014-07-28補(bǔ)充更新

谷延超(1989-)，男，碩士研究生.

P237

：A

：1006-7949(2014)10-0023-04