基于地形的LiDAR數據缺失區域填充方法研究

谷延超，范東明

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

基于地形的LiDAR數據缺失區域填充方法研究

谷延超，范東明

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

針對機載激光點云數據柵格化過程中出現的空白區域進行研究，在對鄰近填充和最低點填充分析的基礎上，提出基于地形的填充方法。對缺失數據邊界進行一維形態學濾波得到其邊界各點地形高，利用邊界點構建不規則三角網進行地形內插填充。實驗結果表明，地形填充法可適應不同原因造成的區域數據缺失，能夠有效保證高程的連續性并可提高濾波精度。

機載LiDAR；柵格化；數據填充；不規則離散點；形態學濾波

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

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

1 點云柵格化

LiDAR獲取的數據為一系列離散的空間三維坐標及反射強度，點云的位置和間隔在空間上都表現為離散性和不規則性。基于規則格網的形態學濾波可直接利用現有的形態學濾波方法進行大數據的處理[1]，故需要對原始的不規則離散數據進行柵格化。

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

2 空白柵格填充

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

當激光腳點的點間距大于影像的空間分辨率時，會產生Ⅰ類數據缺失[2]，可利用限制搜索半徑的鄰近激光點對高程影像和掩膜影像進行填充，搜索半徑可根據設備的掃描方式以及柵格大小來設定，通常可選取1～2倍的柵格大小。

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

2.1 鄰近填充

鄰近填充利用相鄰的地物性質一致進行填充，假設柵格高程影像為I，點云數據集為P，若數據缺失區域為L，則對L內的所有單元格逐一進行最近點搜索，步驟如下：

1)以空白單元格lij為中心，統計搜索半徑r內的所有激光腳點集N。

2)若N為空，則擴大搜索半徑r，重復1)直至腳點集N為非空為止。

3)計算lij與鄰近激光腳點集N中所有腳點的距離，距離最小值所對應的高程即為待填充點的高程。

2.2 最低點填充

文獻[5]中認為：水面往往接近其周邊范圍內的最低點，故提出了最低點填充。對消除Ⅰ類數據缺失后的掩膜影像進行空白區域探測，統計所有的空白區域，依次對Ⅱ類數據缺失利用邊界最低點進行填充，步驟如下[5]：

1)對掩膜影像進行形態學腐蝕運算得到腐蝕后的掩膜影像，將原掩膜影像與腐蝕后的掩膜影像相減即可得到數據缺失的邊界。

2)提取 Ⅱ 類數據缺失的邊界中對應高程影像I的最小值，賦給高程影像中對應的 Ⅱ 類數據缺失區域。

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

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

圖1 邊界提取示意圖

通常情況下，針對由航帶無重疊造成的Ⅱ類數據缺失進行鄰近填補為較優，針對水體吸收造成的Ⅱ類數據缺失進行最低點填補較優[5]，然而在現實數據處理過程中，無法通過空白區域判斷其產生原因，故不能按照上述分類進行填充。

2.3 地形填充

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

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

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

針對 Ⅱ 類數據缺失利用區域邊界，通過一維形態學濾波得到其邊界各點地形，然后利用邊界地形基于TIN對空白區域進行填充。地形填充既能統一不同原因引起的 Ⅱ 類數據缺失，又能保證填充區域內部和外部地形的連續性。地形填充的具體流程如下：

1)對掩膜影像進行形態學腐蝕運算得到腐蝕后的掩膜影像，將原掩膜影像與腐蝕后的掩膜影像相減即可得到數據缺失的邊界。

2)對提取的邊界進行編碼連接。由于提取的邊界點無拓撲關系，而一維形態學濾波是根據周邊地形確定各點高程，故將相鄰的邊界點進行編碼連接構成一維向量，使得相鄰的邊界點在一維向量中同樣緊鄰。

3)邊界點一維形態學濾波。根據實驗區域的先驗知識，確定最大建筑物的尺寸。以最大建筑物的尺寸為濾波結構元的大小，依次得到各邊界點的地形。圖2所示為某一數據缺失區域提取的邊界地形，一維形態學濾波較好地提取了邊界各點的地形。

圖2 某數據缺失區域邊界地形提取示意圖

4)邊界點構TIN并進行內插填充。利用邊界點的像素坐標構建不規則三角形，對空白區域的所有柵格進行判斷其所屬的三角形，然后利用不規則三角形頂點的地形高進行內插求得各空白柵格的地形。

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

3 評價分析

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

3.1 目視分析

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

圖3 數據缺失區域及其填充影像

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

3.2 定量分析

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

濾波質量通過三類誤差進行定量評價，其中第1類誤差(Type Ⅰ)是地面點被誤分為地物點的百分比；第2類誤差(TypeⅡ)是地物點被誤分為地面點的百分比；第3類誤差為總誤差，即被錯分的點占整個數據的百分比[8-9]。

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

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

表1 樣本的三類誤差統計 %

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

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

4 結束語

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

[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]沈晶，劉紀平，林祥國．用形態學重建方法進行機載LiDAR 數據濾波[J]．武漢大學學報:信息科學版，2011，36 (2)：167-176.

[3]李峰，崔希民，袁德寶，等．改進坡度的LiDAR 點云形態學濾波算法[J]．大地測量與地球動力學，2012，32 (5)：128-132.

[4]周曉明，馬秋禾，許嘵亮，等．LIDAR點云濾波算法分析—以ISPRS測試實驗為參考 [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]羅伊萍，姜挺，王鑫,等．基于數學形態學的LiDAR 數據濾波新方法[J]．測繪通報，2011(3)：15-19.

[7]隋立春，張熠斌，柳艷，等．基于改進的數學形態學算法的LiDAR點云數據濾波[J]．測繪學報，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]程曉光，黃先鋒，張帆.機載LiDAR數據的城區樹木點提取方法[J].測繪科學，2014，39(3)：52-56.

[責任編輯：劉文霞]

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補充更新

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

P237

：A

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