王 刃，朱新慧，江振治

(1. 信息工程大學(xué) 導(dǎo)航與空天目標(biāo)學(xué)院，河南 鄭州 450052; 2. 西安測(cè)繪研究所，陜西 西安 710054)

一、引 言

機(jī)載LiDAR數(shù)據(jù)濾波技術(shù)是目前的一個(gè)研究熱點(diǎn)與難點(diǎn)?，F(xiàn)有的機(jī)載LiDAR數(shù)據(jù)濾波技術(shù)大致可以分為兩類(lèi)：一類(lèi)是數(shù)據(jù)驅(qū)動(dòng)式濾波算法，另一類(lèi)是模型驅(qū)動(dòng)式濾波算法。數(shù)據(jù)驅(qū)動(dòng)式濾波算法主要有基于數(shù)學(xué)形態(tài)學(xué)的濾波算法、基于坡度變化的濾波算法、基于TIN結(jié)構(gòu)的分層濾波算法、逐行雙向標(biāo)識(shí)法(one-dimensional and bidirectional labelling，OBL)濾波算法等。模型驅(qū)動(dòng)式濾波算法主要有基于分層穩(wěn)健線性估計(jì)的濾波算法、基于表面自適應(yīng)的濾波算法、基于樣條插值和區(qū)域增長(zhǎng)的濾波算法等[1-6]。這兩類(lèi)算法在分類(lèi)點(diǎn)云時(shí)，絕大多數(shù)僅僅采用了高程差或坡度作為判斷的依據(jù)，從而只能在地形場(chǎng)景單一，地物表面與地形表面高程差異明顯的區(qū)域，獲得較為理想的濾波結(jié)果。但對(duì)于地形場(chǎng)景復(fù)雜，地形存在突變或地物表面與地形表面高程差異不大的區(qū)域，則很難獲得理想的濾波結(jié)果[7-9]。

針對(duì)目前主要的濾波算子所存在的缺陷，同時(shí)兼顧到機(jī)載LIDAR系統(tǒng)獲取數(shù)據(jù)的基本原理，筆者提出一種濾波算子——改進(jìn)的逐行雙向標(biāo)識(shí)法(progressire OBL，POBL)。該濾波算子可有效克服目前濾波算子的缺點(diǎn)，獲得更高準(zhǔn)確率的分類(lèi)結(jié)果。

二、POBL算法的基本原理

POBL濾波算子是對(duì)OBL算子的改進(jìn)，OBL算法是Shan和Sampath于2005年提出的[10]，該算法在過(guò)濾點(diǎn)云數(shù)據(jù)時(shí)，以剖面為處理單位，并在濾波時(shí)采用坡度與高程值兩個(gè)判斷尺度，以克服單一尺度的不足。但由于坡度和與高程差有很大的相關(guān)性，在某種程度上只反映了對(duì)象的同一屬性，因此該濾波算子只適合連續(xù)變化的地形表面數(shù)據(jù)，對(duì)于復(fù)雜場(chǎng)景，該算子濾波效果仍將無(wú)法得到理想的結(jié)果。

POBL濾波算子針對(duì)OBL濾波算子的不足，在濾波時(shí)同時(shí)采用坡度差值、一維長(zhǎng)度尺寸及高程差作為判斷標(biāo)準(zhǔn)，把剖面分割成具有不同一維長(zhǎng)度的特征，當(dāng)特征長(zhǎng)度大于閾值l2時(shí)，則該特征為地形表面；當(dāng)特征長(zhǎng)度位于區(qū)間[l1，l2]內(nèi)，則為建筑表面；當(dāng)特征長(zhǎng)度小于l1時(shí)，則對(duì)應(yīng)的特征為植被表面(如圖1所示)。

圖1 特征長(zhǎng)度分割區(qū)

POBL濾波算法的具體實(shí)現(xiàn)過(guò)程如下：

1) 尋找每一行中高程值最低的點(diǎn)，把該點(diǎn)標(biāo)為地面點(diǎn)。

2) 從地面點(diǎn)開(kāi)始向兩端進(jìn)行標(biāo)識(shí)，一個(gè)稱為前向標(biāo)識(shí)， 另一個(gè)為后向標(biāo)識(shí)。

3) 為了便于闡述具體的標(biāo)識(shí)過(guò)程，首先解釋幾個(gè)標(biāo)簽過(guò)程中涉及的概念。

a. 點(diǎn)的屬性，即點(diǎn)所在表面的對(duì)象類(lèi)型，包含3種值，即地面點(diǎn)、非地面點(diǎn)，以及未知點(diǎn)。

b. 當(dāng)前點(diǎn)的坡度值及坡度差值。A點(diǎn)的坡度值及坡度差值定義為

(1)

式中，ZA為A點(diǎn)處的高程值；ZB為B點(diǎn)處的高程值；SAB為A、B兩點(diǎn)的歐氏距離。

前向標(biāo)簽過(guò)程：

a. 計(jì)算當(dāng)前點(diǎn)的坡度值并和坡度閾值進(jìn)行比較，如果小于坡度閾值，則根據(jù)前一個(gè)點(diǎn)的屬性，判斷當(dāng)前點(diǎn)的屬性，具體來(lái)說(shuō)分為以下幾種情況：

① 前一點(diǎn)為地面點(diǎn)時(shí)，則當(dāng)前點(diǎn)也為地面點(diǎn)。

② 前一個(gè)點(diǎn)為非地面點(diǎn)時(shí)，則需要對(duì)當(dāng)前點(diǎn)的坡度值的正負(fù)進(jìn)行判斷。當(dāng)坡度值為正時(shí)，則當(dāng)前點(diǎn)和前一個(gè)點(diǎn)位于同一平面上；反之，搜索前面已標(biāo)簽過(guò)的點(diǎn)中離當(dāng)前點(diǎn)最近的地面點(diǎn)，利用該點(diǎn)的高程值，以及坡度值對(duì)當(dāng)前點(diǎn)的高程值進(jìn)行預(yù)測(cè)，計(jì)算預(yù)測(cè)值和實(shí)際高程值之間的差值。當(dāng)差值小于事先設(shè)定的高程差閾值時(shí)，則當(dāng)前點(diǎn)為地面點(diǎn)；否則計(jì)算當(dāng)前點(diǎn)的坡度差值，并和坡度差閾值進(jìn)行比較。當(dāng)小于等于坡度差閾值時(shí)，當(dāng)前點(diǎn)和前一個(gè)點(diǎn)位于同一表面上；反之則當(dāng)前點(diǎn)為未知點(diǎn)。

③ 前一點(diǎn)為未知點(diǎn)時(shí)，判斷當(dāng)前點(diǎn)坡度值的正負(fù)，根據(jù)不同符號(hào)值，做不同的分類(lèi)。

當(dāng)符號(hào)值為負(fù)時(shí)，搜索前面已標(biāo)簽過(guò)的點(diǎn)中離當(dāng)前點(diǎn)最近的地面點(diǎn)，利用該地面點(diǎn)的高程值及坡度值對(duì)當(dāng)前點(diǎn)的高程值進(jìn)行預(yù)測(cè)，并計(jì)算當(dāng)前點(diǎn)實(shí)際高程值和該點(diǎn)處預(yù)測(cè)高程值之間的差，如果差值小于高程差閾值，則當(dāng)前點(diǎn)為地面點(diǎn)。同時(shí)，統(tǒng)計(jì)當(dāng)前點(diǎn)之前的所有未知點(diǎn)的個(gè)數(shù)，根據(jù)未知點(diǎn)的總數(shù)，判斷這些點(diǎn)的屬性：當(dāng)大于長(zhǎng)度閾值l2時(shí)，則把它們分為地面點(diǎn)；反之，則為非地面點(diǎn)。如果高程差值大于高程差閾值，計(jì)算當(dāng)前點(diǎn)的坡度差值，并和坡度差閾值進(jìn)行比較，坡度差值小于坡度差閾值，則當(dāng)前點(diǎn)和前一點(diǎn)位于同一表面上，且為未知點(diǎn)；反之則當(dāng)前點(diǎn)為未知點(diǎn)，統(tǒng)計(jì)當(dāng)前點(diǎn)之前的所有未知點(diǎn)的總數(shù)，并和長(zhǎng)度閾值進(jìn)行比較，如果未知點(diǎn)的總數(shù)大于長(zhǎng)度閾值l2時(shí)，則這些未知點(diǎn)為地面點(diǎn)，反之為非地面點(diǎn)。當(dāng)符號(hào)值為正時(shí)，計(jì)算當(dāng)前點(diǎn)的坡度差值，并和坡度差閾值進(jìn)行比較，小于坡度差值閾值，則當(dāng)前點(diǎn)和前一個(gè)點(diǎn)位于同一表面上，且都為未知點(diǎn)；反之，統(tǒng)計(jì)當(dāng)前點(diǎn)之前的所有未知點(diǎn)的個(gè)數(shù)，并和長(zhǎng)度閾值進(jìn)行比較，如果大于長(zhǎng)度閾值l1，則當(dāng)前點(diǎn)為未知點(diǎn)，該點(diǎn)之前的所有未知點(diǎn)為地面點(diǎn)，小于長(zhǎng)度閾值l1則當(dāng)前點(diǎn)為未知點(diǎn)，該點(diǎn)之前的所有未知點(diǎn)為非地面點(diǎn)。

b. 當(dāng)前點(diǎn)的坡度值大于坡度閾值，則根據(jù)前一點(diǎn)的屬性進(jìn)行判斷，可分為以下幾種情況：

① 前一個(gè)點(diǎn)為地面點(diǎn)，則當(dāng)前點(diǎn)為非地面點(diǎn)。

② 前一個(gè)點(diǎn)為非地面點(diǎn)，計(jì)算當(dāng)前點(diǎn)的坡度差值，并和坡度差閾值進(jìn)行比較，小于坡度差閾值，則當(dāng)前點(diǎn)和前一個(gè)點(diǎn)位于同一表面上；反之則當(dāng)前點(diǎn)為未知點(diǎn)。

③ 前一個(gè)點(diǎn)為未知點(diǎn)時(shí)，計(jì)算當(dāng)前點(diǎn)的坡度差值，并和坡度差閾值進(jìn)行比較，小于坡度差閾值，則當(dāng)前點(diǎn)和前一個(gè)點(diǎn)位于同一平面上，且為未知點(diǎn)；反之則當(dāng)前點(diǎn)為未知點(diǎn)，統(tǒng)計(jì)當(dāng)前點(diǎn)之前的所有未知點(diǎn)的總數(shù)，并和長(zhǎng)度閾值進(jìn)行比較，大于長(zhǎng)度閾值，則當(dāng)前點(diǎn)之前的所有未知點(diǎn)為地面點(diǎn)，否則當(dāng)前點(diǎn)的所有未知點(diǎn)為非地面點(diǎn)。

4) 后向標(biāo)識(shí)過(guò)程和前向標(biāo)簽過(guò)程是一致的。

5) 當(dāng)標(biāo)識(shí)過(guò)程遇到一行的端點(diǎn)時(shí)，如果所有點(diǎn)都被賦予了具體的屬性值時(shí)，則標(biāo)識(shí)過(guò)程結(jié)束；如果還有一些點(diǎn)沒(méi)有被賦予具體的屬性值，則統(tǒng)計(jì)所有剩下的未知點(diǎn)的個(gè)數(shù)，并和長(zhǎng)度閾值進(jìn)行比較，根據(jù)比較結(jié)果賦予相應(yīng)的屬性閾值。

6) 在標(biāo)識(shí)完之后，在剖面中進(jìn)行線性擬合，進(jìn)一步剔除非地面點(diǎn)。

三、試驗(yàn)結(jié)果與分析

試驗(yàn)區(qū)域位于地面平坦的郊區(qū)，主要地物為建筑和植被，其中植被包括樹(shù)木、草坪等，并且草坪表面與地形表面較為貼近，建筑中有部分表面積非常平坦。整塊數(shù)據(jù)包含了1 227 593個(gè)離散點(diǎn)，其中，Z值最大為143.24 m，最小為92.99 m，數(shù)據(jù)采樣密度為1.2點(diǎn)/m2(如圖2—圖3所示)。

圖2 原始點(diǎn)云數(shù)據(jù)

試驗(yàn)首先利用基于小波分層原理的濾波算子、基于漸進(jìn)窗口尺寸的數(shù)學(xué)形態(tài)學(xué)濾波算子，以及分層穩(wěn)健線性估計(jì)濾波法，分別對(duì)試驗(yàn)數(shù)據(jù)進(jìn)行濾波處理，并獲取相應(yīng)的濾波結(jié)果；然后把經(jīng)典濾波算子獲取的濾波結(jié)果和POBL濾波算子獲取的濾波結(jié)果進(jìn)行比較。與此同時(shí)，還利用TerraScan軟件通過(guò)人工編輯的方式獲取了對(duì)應(yīng)區(qū)域的DEM作為參考，并通過(guò)比較其他濾波方式獲取的DEM與參考值之間的差異，獲取對(duì)應(yīng)的統(tǒng)計(jì)結(jié)果，這里的DEM高程差異閾值設(shè)為15 cm，試驗(yàn)結(jié)果如圖4—圖7所示，表1為對(duì)應(yīng)的精度統(tǒng)計(jì)結(jié)果。

圖3 對(duì)應(yīng)的航空照片

圖4 POBL濾波獲取的地面點(diǎn)

圖5 漸進(jìn)窗口尺寸數(shù)學(xué)形態(tài)學(xué)濾波算子獲取的地面點(diǎn)

圖6 小波分層濾波獲取的地面點(diǎn)

圖7 分層穩(wěn)健線性估計(jì)濾波獲取的地面點(diǎn)

試驗(yàn)中各濾波算子的參數(shù)設(shè)置為：POBL算子坡度差閾值為8o，高程差閾值為0.3 m，坡度閾值為10°，長(zhǎng)度閾值為100 m；小波分層濾波算子最大窗口尺寸為100 m，坡度閾值為10°；漸進(jìn)窗口尺寸的數(shù)學(xué)形態(tài)濾波算子高程差值閾值為0.3 m，最大窗口尺寸為100 m，坡度閾值為10°；分層穩(wěn)健線性估計(jì)濾波算法C(0)為1，c為60，坡度閾值為10°。分析對(duì)比圖4—圖7可以發(fā)現(xiàn)，經(jīng)典濾波算法對(duì)試驗(yàn)數(shù)據(jù)進(jìn)行濾波時(shí)，位于幾個(gè)建筑物屋頂表面面積較大的激光腳點(diǎn)出現(xiàn)了錯(cuò)誤分類(lèi)的情況(深灰色圓圈內(nèi))，而POBL濾波算子則有效地避免了這一情況的出現(xiàn)。這主要是因?yàn)榍?種濾波算子受限于窗口尺寸與地形細(xì)節(jié)保留之間的矛盾，當(dāng)濾波窗口尺寸較小時(shí)，就會(huì)導(dǎo)致錯(cuò)誤分類(lèi)；而當(dāng)濾波窗口尺寸較大時(shí)，則會(huì)導(dǎo)致地形細(xì)節(jié)被忽略。而POBL濾波算子則可以同時(shí)兼顧地形細(xì)節(jié)保留與大面積建筑物表面點(diǎn)的有效剔除，各類(lèi)濾波算子的性能統(tǒng)計(jì)見(jiàn)表1。

表1 各類(lèi)濾波算子的性能統(tǒng)計(jì)表

由表1可以發(fā)現(xiàn)，POBL算法獲取的濾波結(jié)果內(nèi)插的DEM精度為0.61 m，而其他3種濾波方法獲取的DEM精度在0.91 m左右。

四、結(jié)束語(yǔ)

POBL濾波算法從原理上基本能夠適應(yīng)各類(lèi)地形背景下的數(shù)據(jù)濾波，特別是對(duì)低矮植被點(diǎn)及地形突變的環(huán)境可以取得較好的濾波結(jié)果。同時(shí)，該算法有效擺脫了對(duì)最大窗口尺寸的依賴，可同時(shí)兼顧地形細(xì)節(jié)保留與大面積建筑表面點(diǎn)的有效剔除。

參考文獻(xiàn)：

[1] 黃先鋒，李卉，王瀟，等.機(jī)載LiDAR數(shù)據(jù)濾波方法評(píng)述[J].測(cè)繪學(xué)報(bào)，2009，38(5)：466-469.

[2] 王金亮，陳聯(lián)君. 激光雷達(dá)點(diǎn)云數(shù)據(jù)的濾波算法評(píng)述[J].遙感技術(shù)與應(yīng)用，2010(25)：632-637.

[3] 王明華，張小紅，曾濤，等. 機(jī)載LiDAR數(shù)據(jù)濾波預(yù)處理方法研究[J].武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2010，35(2)：224-227.

[4] 蔣晶鈺，張祖勛，明英. 復(fù)雜城市環(huán)境的機(jī)載LiDAR點(diǎn)云綠波[J].武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2007，32(5)：402-404.

[5] 成曉倩，趙紅強(qiáng). 基于區(qū)域生長(zhǎng)的LIDAR點(diǎn)云數(shù)據(jù)濾波[J].國(guó)土資源遙感，2008(4)：6-8.

[6] 李卉，李德仁，黃先鋒，等.一種漸進(jìn)加密三角網(wǎng)LiDAR點(diǎn)云濾波的改進(jìn)算法[J].測(cè)繪科學(xué)，2009(34)：39-41.

[7] 張玉方，程新方，歐陽(yáng)平，等. 機(jī)載LiDAR數(shù)據(jù)處理及其應(yīng)用綜述[J].工程地球物理學(xué)報(bào)，2008(5)：119-123.

[8] 熊娜，程新文. LiDAR數(shù)據(jù)特點(diǎn)及其濾波分類(lèi)方法[J].中國(guó)水運(yùn)，2008(8)：151-154.

[9] 刁鑫鵬，吳侃. 改進(jìn)的移動(dòng)窗口曲面擬合法點(diǎn)云數(shù)據(jù)濾波處理[J]. 現(xiàn)代礦業(yè)，2011(6)：59-61.

[10] Jie Shan, Aparajithan Sampath. Urban DEM Generation from Raw LiDAR Data: A Labeling Algorithm and Its Performance[J]. Photogrammetric Engineering & Remote Sensing，2005,71(2)：217-226.