林祥國，張繼賢，寧曉剛，段敏燕，臧 藝

中國測繪科學(xué)研究院，北京 100830

?

融合點(diǎn)、對象、關(guān)鍵點(diǎn)等3種基元的點(diǎn)云濾波方法

林祥國，張繼賢，寧曉剛，段敏燕，臧 藝

中國測繪科學(xué)研究院，北京 100830

基元是影響點(diǎn)云濾波精度和效率的關(guān)鍵因素之一。本文提出了一種基于多基元的三角網(wǎng)漸進(jìn)加密(MPTPD)濾波方法。它包括點(diǎn)云分割、對象關(guān)鍵點(diǎn)提取、基于關(guān)鍵點(diǎn)的對象類別判別3個主要階段，且3個階段的基元分別為點(diǎn)、對象、關(guān)鍵點(diǎn)。使用了4景機(jī)載激光雷達(dá)和攝影測量點(diǎn)云數(shù)據(jù)對MPTPD、三角網(wǎng)漸進(jìn)加密(TPD)、基于對象的三角網(wǎng)漸進(jìn)加密(OTPD)3種濾波方法進(jìn)行了性能測試。試驗(yàn)表明，MPTPD方法具有整體上最優(yōu)的性能：在精度方面，MPTPD與OTPD兩種方法的精度相當(dāng)，MPTPD方法的一類誤差I(lǐng)、總誤差T比TPD的相應(yīng)誤差分別低約22.07%和8.44%；在效率方面，多數(shù)情況下TPD、MPTPD、OTPD方法的效率依次降低，且MPTPD的平均耗時是OTPD平均耗時的57.93%。

濾波；激光雷達(dá)點(diǎn)云；攝影測量點(diǎn)云；對象；三角網(wǎng)

隨著激光雷達(dá)(light detection and ranging，LiDAR)[1]測量、多視影像密集匹配[2]技術(shù)的完善和行業(yè)應(yīng)用的深入，點(diǎn)云濾波的重要性日益突出。本文的點(diǎn)云涉及機(jī)載LiDAR點(diǎn)云和航空、航天多視立體影像密集匹配的點(diǎn)云等3種類型。在點(diǎn)云處理和信息提取領(lǐng)域，濾波是指區(qū)分點(diǎn)云中的地面點(diǎn)和非地面點(diǎn)的過程[1,3]，它是生成數(shù)字高程模型(DEM)、分類、目標(biāo)識別和三維重建的基礎(chǔ)和必經(jīng)的步驟[4]。文獻(xiàn)[1,4-5]對目前眾多點(diǎn)云濾波方法進(jìn)行了系統(tǒng)的介紹。其中，有代表性的方法有三角網(wǎng)(triangular irregular network，TIN)漸進(jìn)加密(TIN progressive densification，TPD)[6-7]、分層穩(wěn)健線性內(nèi)插[8]、坡度濾波[9]、數(shù)學(xué)形態(tài)學(xué)濾波[10]、基于聚類/對象的濾波[11-12]等。已有方法中涉及的基元(基本處理單元)有點(diǎn)[6-7]、對象[11-12]、體素[13]或剖面[14]等多種類型；且后3種基元具有一定的共性，本質(zhì)上是點(diǎn)基元的一種集合和再組織方式，本文僅關(guān)注其中的對象。由于點(diǎn)易受粗差、地形斷裂的負(fù)面影響，而對象比點(diǎn)更能增強(qiáng)點(diǎn)云處理效果[15]，因此基于對象的點(diǎn)云濾波方法[16-21]是研究的一個熱點(diǎn)。然而，與基于點(diǎn)的濾波方法相比，盡管基于對象的濾波方法可以在一定程度上提高濾波精度，但是也存在效率低下的問題[21]。

文獻(xiàn)[22]提出基于多實(shí)體的點(diǎn)云分類方法，在分類的不同階段使用不同的實(shí)體以實(shí)現(xiàn)更優(yōu)的分類效果。借鑒上述策略，本文設(shè)計一種既能繼承基于對象方法的優(yōu)勢、又不顯著降低基于點(diǎn)方法的效率的濾波技術(shù)，即同時提高基于點(diǎn)的TPD方法[6]精度和基于對象的三角網(wǎng)漸進(jìn)加密(object-based TPD，OTPD)方法[21]效率，本文稱之為基于多基元的三角網(wǎng)漸進(jìn)加密(multiple-primitives-based TPD，MPTPD)方法。它有3個創(chuàng)新點(diǎn)：① 使用多基元、而非單一的基元參與運(yùn)算，其中多基元包括點(diǎn)、對象、關(guān)鍵點(diǎn)等3種類型，且在不同階段使用不同類型的基元；② 使用關(guān)鍵點(diǎn)代替對象參與判別，即在核心判別步驟中，使用對象的關(guān)鍵點(diǎn)替代對象進(jìn)行運(yùn)算以提高效率；③ 提出一種簡單、快捷的關(guān)鍵點(diǎn)檢測算法。

特別指出，本文的一個“對象”指“點(diǎn)云分割后具有同一標(biāo)號的點(diǎn)集”，“關(guān)鍵點(diǎn)”又是對象點(diǎn)集的一個子集，即關(guān)鍵點(diǎn)本質(zhì)上仍然是原始點(diǎn)云中的點(diǎn)，而非額外創(chuàng)造的，但是關(guān)鍵點(diǎn)具有特殊性。另外，處理一個“對象”，可以通過處理該對象包含的點(diǎn)集來實(shí)現(xiàn)，也可以通過處理“關(guān)鍵點(diǎn)”來實(shí)現(xiàn)。

1 本文的濾波方法

MPTPD方法包括基于表面生長的點(diǎn)云分割、對象關(guān)鍵點(diǎn)提取、基于關(guān)鍵點(diǎn)的對象類別判別等3個主要步驟。整體技術(shù)框架如圖1所示；圖2展示了某一點(diǎn)云各個處理步驟的效果，文中數(shù)字“1”代表“非地面點(diǎn)類”，數(shù)字“2”代表“地面點(diǎn)類”。

圖1 本文方法的整體技術(shù)流程圖Fig.1 The whole work flow of the proposed method

1.1 基于表面生長的點(diǎn)云分割

點(diǎn)云分割是對點(diǎn)云數(shù)據(jù)中每個點(diǎn)按照一定的判別規(guī)則進(jìn)行標(biāo)號的過程。分割后，滿足同一規(guī)則的點(diǎn)集被賦予同一標(biāo)號，且每一點(diǎn)集稱為一個對象。本文判別的是3D空間中鄰近且共平面的點(diǎn)。另外，不滿足上述判別規(guī)則的孤立點(diǎn)、鄰近點(diǎn)數(shù)目不足的點(diǎn)、共平面性差的點(diǎn)亦會被標(biāo)號。鑒于表面生長[23]算法具有所需參數(shù)少、分割效果好、普適性好的特點(diǎn)，本文使用它對點(diǎn)云進(jìn)行分割，其主要步驟如下。

第1步，估計法向量和殘差，處理過程如下。

(1) 加載點(diǎn)云數(shù)據(jù)，并將所有點(diǎn)的類別標(biāo)記為“1”、標(biāo)號狀態(tài)記為“未分割”，設(shè)共計有n個點(diǎn)。

(2) 建立點(diǎn)云的三維kd-tree[24]空間索引。

(3) 逐一處理每一個點(diǎn)，即對第i(i=0,2,…，n-1)個點(diǎn)，首先利用kd-tree求取其k個最臨近點(diǎn)，然后利用特征值法[25]求當(dāng)前點(diǎn)及k個鄰近點(diǎn)構(gòu)成點(diǎn)集的擬合平面方程，即可確定第i個點(diǎn)的法向量φi及其殘差λi。

第2步，進(jìn)行區(qū)域生長。

圖2 表面生長過程的示意圖Fig.2 Illustration of the surface growing process

(1) 記“區(qū)域標(biāo)記號”從0開始。

(2) 檢查“未分割”點(diǎn)集中點(diǎn)的數(shù)量，如果數(shù)量為0，則轉(zhuǎn)到步驟(6)；否則，接著從“未被分割”的點(diǎn)集中，尋找出殘差λ最小的點(diǎn)，以該點(diǎn)為種子點(diǎn)并將該點(diǎn)壓入一個種子點(diǎn)的隊(duì)列，且將該點(diǎn)的處理狀態(tài)標(biāo)記為“未處理”，開始進(jìn)行區(qū)域生長。

(3) 取種子點(diǎn)隊(duì)列中第一個“未處理”的種子點(diǎn)，利用kd-tree求取該種子點(diǎn)的k個最臨近點(diǎn)。

(4) 逐一對于每一個臨近點(diǎn)進(jìn)行判別。如果臨近點(diǎn)已經(jīng)被賦予分割號，則不予以處理；另外，若臨近點(diǎn)已經(jīng)在種子點(diǎn)的隊(duì)列中，則不予以處理；接著，分別按照法向量間角度差異和距離差異的規(guī)則進(jìn)行當(dāng)前種子點(diǎn)和該鄰近點(diǎn)的相似性判別。如果該鄰近點(diǎn)同時滿足兩個相似性的條件，則將該鄰近點(diǎn)加入到種子點(diǎn)隊(duì)列中；反之，如果該鄰近點(diǎn)沒有同時滿足兩個相似性的條件，則不予以處理。鄰近點(diǎn)判別完畢后，將該種子點(diǎn)的處理狀態(tài)標(biāo)記為“已處理”。

(5) 順序檢查種子點(diǎn)隊(duì)列中是否有“未處理”的點(diǎn)。如果有，返回步驟(3)；否則，將種子點(diǎn)隊(duì)列中的點(diǎn)集的標(biāo)號記為“區(qū)域標(biāo)記號”，狀態(tài)記為“已分割”，同時“區(qū)域標(biāo)記號”自增1，清空種子點(diǎn)隊(duì)列后返回步驟(2)。

(6) 結(jié)束。

經(jīng)過上述分割后，任意一點(diǎn)被劃到一個對象，但是部分對象的點(diǎn)的數(shù)量較少。個別情況下，一個對象僅僅包含一個點(diǎn)。圖3(a)展示了某機(jī)載LiDAR點(diǎn)云數(shù)據(jù)，圖3(b)展示了其分割效果，其中地面點(diǎn)被聚為若干個對象，多數(shù)地面對象包含點(diǎn)的數(shù)量較多；一個建筑物可以被聚為一個或者若干個對象，這與其類型、點(diǎn)云密度、精度等多個因素有關(guān)；孤立的植被點(diǎn)、粗差點(diǎn)也往往被分割為一個對象。

圖3 本文提出濾波方法的關(guān)鍵步驟處理效果示意圖Fig.3 Illustration of the process of the proposed filtering method

1.2 對象關(guān)鍵點(diǎn)的提取

本文的關(guān)鍵點(diǎn)包括外輪廓點(diǎn)、內(nèi)特征點(diǎn)、最高點(diǎn)和最低點(diǎn)。

第1步，計算每個對象包含點(diǎn)的數(shù)量。如果數(shù)量不大于經(jīng)驗(yàn)閾值4，則將該對象的點(diǎn)集記為對象的關(guān)鍵點(diǎn)。反之，進(jìn)入第2步。

第2步，識別每個對象的關(guān)鍵點(diǎn)。圖4展示了某對象關(guān)鍵點(diǎn)檢測的主要過程。提取的基本原理是僅利用某一對象點(diǎn)集的水平坐標(biāo)信息生成TIN。該TIN中，處于邊緣的三角形僅有兩個三角形通過邊相鄰；而處于非邊緣的三角形有3個三角形通過邊相鄰。因此，可以通過一個三角形的一邊為鄰邊的鄰接三角形的數(shù)量來判斷該三角形是否處于邊緣。處于邊緣的三角形涉及的3個頂點(diǎn)記為“外輪廓點(diǎn)”，如圖4(b)所示。接著，刪除該TIN中的短邊，并以刪除短邊后的TIN為索引進(jìn)行連通區(qū)域分析(connectedcomponentanalysis)。則該對象被分割為若干子對象。若子對象包含的點(diǎn)的數(shù)量大于經(jīng)驗(yàn)閾值4，則該子對象的點(diǎn)集被認(rèn)為是“內(nèi)特征點(diǎn)”，如圖4(c)所示。注意上述關(guān)鍵點(diǎn)無需有序排列，這與Alphashape[26]算法有著顯著的差別。第2步具體內(nèi)容如下。

(1) 建立某一對象的二維TIN。

(2) 檢測“外輪廓點(diǎn)”。通過上述鄰接三角形數(shù)量的規(guī)則識別“外輪廓點(diǎn)”，如圖4(b)所示。

(3) 刪除短邊。假設(shè)輸入原始點(diǎn)云的平均點(diǎn)間距為g(單位：m)，且g已知。刪除TIN中二維邊緣長度小于經(jīng)驗(yàn)閾值3g的邊。

(4) 檢測“內(nèi)特征點(diǎn)”。通過上述連通區(qū)域分析獲取的子對象包含點(diǎn)的數(shù)量的規(guī)則識別“內(nèi)特征點(diǎn)”，如圖4(c)所示。

(5) 檢測最高點(diǎn)、最低點(diǎn)。另外，外輪廓點(diǎn)、內(nèi)特征點(diǎn)、最高點(diǎn)、最低點(diǎn)不可重復(fù)。如有重復(fù)，則只保留其中一個。

圖4(a)展示了某一地面對象包含321 998個點(diǎn)，其外輪廓點(diǎn)、內(nèi)特征點(diǎn)和關(guān)鍵點(diǎn)分別如圖4(b)、(c)和(d)所示，圖4(d)只包含19 875個點(diǎn)。從圖4(a)和(d)中的DEM看，盡管構(gòu)建DEM的點(diǎn)的數(shù)量差別懸殊，但是DEM的表達(dá)效果卻趨于一致。對兩個DEM，不僅最高、最低值一致，且相應(yīng)像素值之差的絕對值的平均值和標(biāo)準(zhǔn)差分別為0.05和0.01，這反映了提取的關(guān)鍵點(diǎn)既能顯著地減少點(diǎn)的數(shù)量、又能逼近真實(shí)的對象原始形態(tài)。另外，圖2(c)展示了圖2(a)中點(diǎn)云的關(guān)鍵點(diǎn)檢測結(jié)果，其中原始點(diǎn)云包含826 416個點(diǎn)、而關(guān)鍵點(diǎn)只包含64 257個點(diǎn)，關(guān)鍵點(diǎn)數(shù)量只占原始點(diǎn)數(shù)量的7.78%。

圖4 對象關(guān)鍵點(diǎn)檢測的示意圖Fig.4 Illustration of detection of the key points of an object

1.3 基于關(guān)鍵點(diǎn)的對象類別判別

經(jīng)典TPD方法的運(yùn)算過程中，TIN構(gòu)建和點(diǎn)類別判別占整個濾波時間的比重很大[27]。本文采用關(guān)鍵點(diǎn)替代對象的目的是能同時顯著地減少參與TIN構(gòu)建的、參與判別的點(diǎn)的數(shù)量以提高效率，又可使構(gòu)建的TIN盡可能地逼近區(qū)域真實(shí)的DEM以確保精度。

本節(jié)是一個迭代過程。共涉及4個參數(shù)：最大建筑物長度b(單位：m)、最大角度閾值θ(單位：°)、最大距離閾值d(單位：m)和最大地形角度閾值t(單位：°)。具體過程如下。

第1步：格網(wǎng)劃分。求點(diǎn)云在XOY平面上的最小外包矩形，并在XOY平面上對該最小外包矩形進(jìn)行格網(wǎng)劃分、且格網(wǎng)的尺寸為b×b。

第2步：將全部對象處理狀態(tài)均標(biāo)記為“未處理”。

第3步：選擇初始地面種子點(diǎn)。逐一選擇每一個格網(wǎng)的地面種子點(diǎn)。即對每一個格網(wǎng)，找到格網(wǎng)中高程值最低的點(diǎn)所在的對象。如果對象的面積小于4.00m2，則繼續(xù)找到高程值次低的點(diǎn)所在的對象直至找到面積大于4.00m2的對象。將該對象的關(guān)鍵點(diǎn)作為該格網(wǎng)的地面種子點(diǎn)，且該對象的類別被標(biāo)記為“2”、處理狀態(tài)被標(biāo)記為“已處理”。

第4步：構(gòu)建初始地面種子點(diǎn)的TIN。該TIN代表該區(qū)域初始的DEM。

第5步：迭代的判別對象類別。子步驟包括：

(1) 迭代次數(shù)記為0。

(2) 以對象為基本處理單元，逐一通過每一個“未處理”對象的關(guān)鍵點(diǎn)的判別，實(shí)現(xiàn)該對象類別的判別。

對每一個“未處理”對象逐一判別其“未處理”關(guān)鍵點(diǎn)，找到該關(guān)鍵點(diǎn)落入的三角形，計算該關(guān)鍵點(diǎn)到三角形構(gòu)成的平面的距離及該關(guān)鍵點(diǎn)到三角形3點(diǎn)的夾角，并找出3個夾角中的最大角。進(jìn)行下述判別：如果同時滿足距離小于d、最大夾角小于θ，則認(rèn)為該關(guān)鍵點(diǎn)是地面點(diǎn)，將該關(guān)鍵點(diǎn)的類別號標(biāo)記為“2”，處理狀態(tài)標(biāo)記為“已處理”，繼續(xù)處理下一個“未處理”關(guān)鍵點(diǎn)；否則，檢查三角形構(gòu)成的平面的傾角，進(jìn)行下述判別：如果傾角小于，繼續(xù)處理下一個“未處理”關(guān)鍵點(diǎn)；反之，將當(dāng)前關(guān)鍵點(diǎn)以所在三角形的最高點(diǎn)為中心做一個鏡像點(diǎn)(參考文獻(xiàn)[6，12])，且該鏡像點(diǎn)的高等于該關(guān)鍵點(diǎn)。對該鏡像點(diǎn)進(jìn)行類似的判別。如果該鏡像點(diǎn)被判別為“2”，則將該關(guān)鍵點(diǎn)的類別號標(biāo)記為“2”，處理狀態(tài)標(biāo)記為“已處理”；繼續(xù)處理下一個“未處理”關(guān)鍵點(diǎn)。

判別完畢，統(tǒng)計該對象的關(guān)鍵點(diǎn)的數(shù)量、關(guān)鍵點(diǎn)屬于地面點(diǎn)的比例。如果該比例大于50%，將該對象的處理狀態(tài)標(biāo)記為“已處理”，類別為“2”；否則，將該對象的處理狀態(tài)重新標(biāo)記為“未處理”，類別為“1”。

(3) 利用新識別的屬于地面點(diǎn)的關(guān)鍵點(diǎn)更新TIN，同時迭代次數(shù)自增1。

(4) 重復(fù)上述步驟(1)至(3)，繼續(xù)執(zhí)行直至迭代次數(shù)達(dá)到經(jīng)驗(yàn)閾值5，或者沒有新識別的地面點(diǎn)則停止迭代。

2 試驗(yàn)與分析

基于VisualStudio2010C++集成開發(fā)環(huán)境實(shí)現(xiàn)了本文提出的MPTPD方法，同時對TPD[6]、OTPD[21]兩種方法進(jìn)行性能比較，上述3種方法均采用串行計算，未采用并行計算技術(shù)。其中，TPD包括低位粗差點(diǎn)剔除和1.3節(jié)描述的5個主要步驟，其基元為點(diǎn)；OTPD包括1.1節(jié)和1.3節(jié)兩個相似部分，其基元為對象。為了增加效率的可比性，盡管3種濾波方法的基元的不同，但相似步驟涉及的算法一致。試驗(yàn)平臺的配置：ThinkPadW520筆記本，CPU為Intel酷睿i7-2760QM2.4GHz，內(nèi)存2.98GB，裝配WindowsXP系統(tǒng)。

2.1 試驗(yàn)數(shù)據(jù)及結(jié)果

本文共使用了4個場景的點(diǎn)云開展試驗(yàn)(圖5)，它們的基本信息見表1。前兩個場景的點(diǎn)云為開放的機(jī)載LiDAR數(shù)據(jù)；后兩個為攝影測量點(diǎn)云，其中，第3個場景的原始影像由TrimbleGermanyGmbH公司免費(fèi)提供、點(diǎn)云由中國測繪科學(xué)研究院的PixelGrid軟件生成，第4個場景的點(diǎn)云由德國宇航局免費(fèi)提供。另外，試驗(yàn)數(shù)據(jù)1為國際攝影測量與遙感協(xié)會第三委員會提供的測試數(shù)據(jù)CSite1，該數(shù)據(jù)位于德國，如圖5(a)所示；試驗(yàn)數(shù)據(jù)2由IEEEGRSLFusionContest2013提供，該數(shù)據(jù)位于美國休斯敦大學(xué)附近，本文截取了原始點(diǎn)云的一部分，如圖5(b)所示；試驗(yàn)數(shù)據(jù)3對應(yīng)的斜影像由天寶公司AOS系統(tǒng)獲取，該數(shù)據(jù)位于德國柏林市和波茨坦市附近，本文截取了由影像生成的攝影測量點(diǎn)云的一部分，如圖5(c)所示；試驗(yàn)數(shù)據(jù)4由GeoEye-1的立體影像對生成，該數(shù)據(jù)位于德國的慕尼黑市，本文截取了該攝影測量點(diǎn)云的一部分，如圖5(d)所示。

對4個試驗(yàn)數(shù)據(jù)進(jìn)行濾波時，使用的相關(guān)參數(shù)的值見表2。其中，TPD、OTPD和MPTPD3種方法對1.3節(jié)的參數(shù)b、θ、d、t采用了相同的參數(shù)值，OTPD和MPTPD兩種方法對1.1節(jié)的參數(shù)k、d、r采用了相同的參數(shù)值。鑒于篇幅的原因未展示3種方法的濾波結(jié)果。后續(xù)試驗(yàn)分析表明3種濾波方法均能正確地區(qū)分多數(shù)的地面點(diǎn)和非地面點(diǎn)，且OTPD和MPTPD兩種方法的濾波效果相當(dāng)、并優(yōu)于TPD方法的濾波效果。

2.2 精度評價

本文采用文獻(xiàn)[4]中的一類誤差I(lǐng)(將地面點(diǎn)錯分為非地面點(diǎn)的數(shù)量占地面點(diǎn)數(shù)量的比例)、二類誤差I(lǐng)I(將非地面點(diǎn)錯分為地面點(diǎn)的數(shù)量占非地面點(diǎn)數(shù)量的比例)和總誤差T(錯分點(diǎn)數(shù)量占全部點(diǎn)數(shù)量的比例)3個指標(biāo)定量衡量濾波精度。同時，使用了人工半自動解譯的方式識別了4個試驗(yàn)數(shù)據(jù)的地面點(diǎn)和非地面點(diǎn)，并將每個試驗(yàn)數(shù)據(jù)的人工識別結(jié)果作為真值計算濾波方法的誤差。4個試驗(yàn)數(shù)據(jù)的3類誤差值見表3。

圖5 4個試驗(yàn)數(shù)據(jù)Fig. 5 The four testing datasets

表3中，TPD方法的Ⅰ、Ⅱ、T3類誤差的平均值分別為27.09%、2.27%和12.52%，OTPD方法相關(guān)誤差的平均值分別為4.96%、2.75%和4.04%，MPTPD方法相關(guān)誤差的平均值分別為5.03%、2.78%和4.08%。數(shù)字說明OTPD方法與MPTPD方法的各類誤差均十分接近。且，一類誤差I(lǐng)和總誤差T均呈現(xiàn)TPD>OTPD≈MPTPD的趨勢，而二類誤差I(lǐng)I呈現(xiàn)TPD

表1 4個試驗(yàn)數(shù)據(jù)的基本信息

表2 4個試驗(yàn)數(shù)據(jù)中3種濾波方法的相關(guān)參數(shù)取值

表3 4個試驗(yàn)數(shù)據(jù)中3種濾波方法3類誤差的統(tǒng)計值

Tab.3 The vales of three types of errors of the three filtering methods for the four testing datasets

試驗(yàn)數(shù)據(jù)誤差類型TPD/(%)OTPD/(%)MPTPD/(%)1Ⅰ34．7411．3511．36Ⅱ2．394．214．16T17．927．647．622Ⅰ22．466．816．98Ⅱ3．020．790．77T13．354．044．073Ⅰ23．711．291．33Ⅱ1．805．745．81T9．404．204．264Ⅰ27．450．390．43Ⅱ1．870．240．36T9．400．280．38

2.3 效率評價

本文采用時間花費(fèi)來衡量濾波效率。為此，分別記錄了3種方法對4個試驗(yàn)數(shù)據(jù)進(jìn)行濾波的各個階段耗時及總耗時，具體的統(tǒng)計數(shù)據(jù)見表4。其中，將濾波過程劃分為5個階段：低位粗差點(diǎn)剔除、基于表面生長的點(diǎn)云分割(1.1節(jié))、對象關(guān)鍵點(diǎn)提取(1.2節(jié))、基于關(guān)鍵點(diǎn)的對象類別判別(1.3節(jié))前3步和后兩步，在表4中分別稱為階段1、階段2、階段3、階段4、階段5。

表4的統(tǒng)計數(shù)據(jù)表明，在4個試驗(yàn)數(shù)據(jù)中，有3個出現(xiàn)了3種濾波方法的總效率均呈現(xiàn)TPD>MPTPD>OTPD的規(guī)律。以試驗(yàn)數(shù)據(jù)3為例，TPD、MPTPD、OTPD 3種方法的總耗時由少到多依次為約299 s、415 s、866 s。但在第1個試驗(yàn)數(shù)據(jù)中，MPTPD的效率高于TPD、OTPD。整體上，TPD效率最高，MPTPD次之，OTPD最低。表4還表明TPD、MPTPD、OTPD 3種方法在每個試驗(yàn)數(shù)據(jù)上的總耗時平均分別為248.76 s、355.26 s、613.30 s。如果以效率最慢的OTPD的基準(zhǔn)，TPD和MPTPD的效率分別是OTPD的2.47倍、1.73倍。

而且，表4還表明每種濾波方法的各個階段的耗時的比例也有著顯著的差別，表現(xiàn)為：

(1) TPD方法中，階段1、階段4和階段5這3個階段的耗時占總耗時比例的平均值分別為84.61%、6.44%、8.95%，可見，階段1和階段5這兩個階段占了TPD方法總耗時的絕大部分比例，其中粗差剔除的相關(guān)比例最大、且點(diǎn)云中的粗差越復(fù)雜相應(yīng)的比例越大。例如，試驗(yàn)數(shù)據(jù)1的粗差多且多樣，粗差剔除的時間花費(fèi)占總時間花費(fèi)的95.19%，是4個試驗(yàn)數(shù)據(jù)中比例最大的，這是在第1個試驗(yàn)中MPTPD的效率高于TPD的原因。

(2) OTPD方法中，階段2、階段4、階段5等3個階段的耗時占總耗時比例的平均值分別為45.62%、10.19%、44.87%。可見，OTPD的3個階段中，除了第2個階段外，其他兩個階段均比較耗時、且占總耗時的比例相當(dāng)。

(3) MPTPD方法中，階段2、階段3、階段4、階段5等4個階段的耗時占總耗時比例的平均值分別為75.47%、20.86%、2.87%、0.80%?？梢姡琈PTPD的4個階段中，前兩個階段的累計耗時占了總耗時的絕大部分，而后兩個階段的累計耗時占總耗時的比例極小。這是MPTPD與OTPD盡管在階段2的耗時相同、但是MPTPD的總耗時顯著的低于OTPD的原因。

另外，表4展示了一個很有興趣的現(xiàn)象。3種濾波方法均有階段4、階段5兩個階段，但由于基元的不同導(dǎo)致兩個階段(尤其是階段5)的效率有著顯著差別。例如，第4個試驗(yàn)數(shù)據(jù)中，TPD和OTPD兩種方法的階段4的耗時分別為9.45 s、109.70 s，但MPTPD的相應(yīng)耗時僅為7.78 s；TPD、OTPD、MPTPD 3種方法在階段5的耗時分別為13.79 s、342.64 s、2.09 s，即MPTPD方法在此階段的耗時僅為TPD的15.16%、OTPD的0.61%。其他3個試驗(yàn)數(shù)據(jù)亦表現(xiàn)出類似的規(guī)律。這證明了基元對濾波效率有著顯著的影響，MPTPD方法的關(guān)鍵點(diǎn)顯著提高了其后續(xù)階段的效率。但是，與TPD方法相比，MPTPD方法的點(diǎn)云分割和關(guān)鍵點(diǎn)提取耗費(fèi)了較多的時間，因此該方法的整體效率低于TPD方法。

2.4 分析與討論

文獻(xiàn)[4]指出TPD方法對具有不同場景復(fù)雜度的點(diǎn)云數(shù)據(jù)均具有較高的濾波精度。本文的統(tǒng)計數(shù)據(jù)表明，TPD方法的總誤差T平均值約12.52%，精度較高，符合既有結(jié)論。另外，表2表明TPD方法所需的4個參數(shù)的取值對場景的變化不是很敏感，且具有顯著的物理意義，根據(jù)實(shí)際情況微調(diào)參數(shù)取值即可。但TPD方法存在對低位粗差、地形斷裂敏感的問題，因此該濾波方法仍然存在一定的誤差。表3表明TPD方法的一類誤差I(lǐng)平均值約27.09%，顯著地高于其他兩種濾波方法的相關(guān)誤差值。

OTPD方法是對TPD方法的改進(jìn)，具有對低位粗差、地形斷裂不敏感的優(yōu)勢，但是比較耗時。統(tǒng)計數(shù)據(jù)表明，與TPD相比，OTPD方法的一類誤差I(lǐng)、總誤差T比TPD的分別低約22.13%、8.48%，但耗時是TPD的2.47倍。

本文提出的MPTPD方法，既有與OTPD方法相當(dāng)?shù)臑V波精度，又有更高的效率。統(tǒng)計數(shù)據(jù)表明，MPTPD方法與OTPD方法的一類誤差Ⅰ、二類誤差Ⅱ、總誤差T的差值絕對值分別為0.07%、0.03%、0.04%，但是MPTPD方法的平均耗時卻僅有OTPD方法的約58%。與OTPD方法的濾波精度相當(dāng)、但效率顯著提升的原因在于，MPTPD方法中的關(guān)鍵點(diǎn)既能逼近原始點(diǎn)云的三維形態(tài)、又能顯著地減少參與后續(xù)判別的計算量。

另外，與TPD方法相比，MPTPD和OTPD兩種方法需要額外的3個參數(shù)。但是，表2表明4個試驗(yàn)數(shù)據(jù)中k和r兩個參數(shù)的取值可相同，而α的取值在3個試驗(yàn)數(shù)據(jù)亦相同。這表明額外的3個參數(shù)亦“具有顯著的物理意義，根據(jù)實(shí)際情況微調(diào)參數(shù)取值即可”。

表4 4個試驗(yàn)數(shù)據(jù)中3種濾波方法的時間花費(fèi)

3 結(jié) 論

點(diǎn)云是一種新型的數(shù)據(jù)源，其數(shù)據(jù)處理方法亟待研究。濾波是點(diǎn)云數(shù)據(jù)處理的一個必要的關(guān)鍵環(huán)節(jié)。目前，多數(shù)濾波方法采用單一的基元，但采用單一基元的濾波方法很難平衡濾波精度和濾波效率。為此，本文提出了MPTPD方法，該方法在濾波的不同階段采用了不同的基元。其涉及的基元包括點(diǎn)、對象、關(guān)鍵點(diǎn)等3種，且濾波的基本原理與TPD、OTPD方法相似。采用4個有代表性的點(diǎn)云數(shù)據(jù)進(jìn)行了試驗(yàn)。試驗(yàn)表明，MPTPD方法具有整體上最優(yōu)的性能。其中，精度方面，MPTPD與OTPD兩種方法的精度相當(dāng)，MPTPD方法的一類誤差I(lǐng)、總誤差T比TPD的相應(yīng)誤差低分別約22.07%、8.44%；效率方面，多數(shù)情況下TPD、MPTPD、OTPD方法的效率依次降低，但少數(shù)情況下MPTPD的效率最高，且MPTPD的平均耗時是OTPD平均耗時的57.93%。筆者下一步的研究圍繞兩個方面開展：①采用并行計算技術(shù)提高濾波效率；②探索多基元的點(diǎn)云分類。

[1] 黃先鋒, 李卉, 王瀟, 等. 機(jī)載LiDAR數(shù)據(jù)濾波方法評述[J]. 測繪學(xué)報, 2009, 38(5): 466-469.

HUANG Xianfeng, LI Hui, WANG Xiao, et al. Filter Algorithms of Airborne LiDAR Data: Review and Prospects[J]. Acta Geodaetica et Cartographica Sinica, 2009, 38(5): 466-469.

[2] 王競雪, 朱慶, 王偉璽. 多匹配基元集成的多視影像密集匹配方法[J]. 測繪學(xué)報, 2013, 42(5): 691-698.

WANG Jingxue, ZHU Qing, WANG Weixi. A Dense Matching Algorithm of Multi-view Image Based on the Integrated Multiple Matching Primitives[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(5): 691-698.

[3] 張小紅. 機(jī)載激光雷達(dá)測量技術(shù)理論與方法[M]. 武漢: 武漢大學(xué)出版社, 2007.

ZHANG Xiaohong. Theories, Methods of Airborne LiDAR Technique[M]. Wuhan: Wuhan University Press, 2007.

[4] 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.

[5] MENG Xuelian, CURRIT N, ZHAO Kaiguang. Ground Filtering Algorithms for Airborne LiDAR Data: A Review of Critical Issues[J]. Remote Sensing, 2010, 2(3): 833-860.

[6] AXELSSON P E. DEM Generation from Laser Scanner Data Using Adaptive TIN Models[J]. International Archives of Photogrammetry and Remote Sensing, 2000, 33: 110-117.

[7] 隋立春, 張熠斌, 張碩, 等. 基于漸進(jìn)三角網(wǎng)的機(jī)載LiDAR點(diǎn)云數(shù)據(jù)濾波[J]. 武漢大學(xué)學(xué)報(信息科學(xué)版), 2011, 36(10): 1159-1163. SUI Lichun, ZHANG Yibin, ZHANG Shuo, et al. Filtering of Airborne LiDAR Point Cloud Data Based on Progressive TIN[J]. Geomatics and Information Science of Wuhan University, 2011, 36(10): 1159-1163.

[8] BRIESE C, PFEIFER N, DORNINGER P. Applications of the Robust Interpolation for DTM Determination[J]. International Archives of Photogrammetry and Remote Sensing, 2002, 34: 55-61.

[9] SITHOLE G. Filtering of Laser Altimetry Data Using a Slope Adaptive Filter[J]. International Archives of Photogrammetry and Remote Sensing, 2001, 34: 203-210.

[10] ZHANG Keqi, CHEN Shuching, WHITMAN D, et al. A Progressive Morphological Filter for Removing Nonground Measurements from Airborne LiDAR Data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2003, 41(4): 872-882.

[12] ZHANG Jixian, LIN Xiangguo. Filtering Airborne LiDAR Data By Embedding Smoothness-constrained Segmentation in Progressive TIN Densification[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2013, 81: 44-59.

[13] 唐菲菲, 劉經(jīng)南, 張小紅, 等. 基于體素的森林地區(qū)機(jī)載LiDAR數(shù)據(jù)DTM提取[J]. 北京林業(yè)大學(xué)學(xué)報, 2009, 31(1): 55-59.

TANG Feifei, LIU Jingnan, ZHANG Xiaohong, et al. A Voxel-based Filtering Algorithm for DTM Data Extraction in Forest Areas[J]. Journal of Beijing Forestry University, 2009, 31(1): 55-59.

[14] 鄭輯濤, 張濤. 基于可變半徑圓環(huán)和B樣條擬合的機(jī)載LiDAR點(diǎn)云濾波[J]. 測繪學(xué)報, 2015, 44(12): 1359-1366. DOI: 10.11947/j.AGCS.2015.20140514.

ZHENG Jitao, ZHANG Tao.Filtering of Airborne LiDAR Point Cloud Based on Variable Radius Circle and B-spline Fitting[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(12): 1359-1366. DOI: 10.11947/j.AGCS.2015.20140514.

[15] FILIN S, PFEIFER N. Segmentation of Airborne Laser Scanning Data Using A Slope Adaptive Neighborhood[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2006, 60(2): 71-80.

[16] LOHMANN P. Segmentation and Filtering of Laser Scanner Digital Surface Models[J]. International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences, 2002, 34(Pt 3): 311-315.

[17] CHEN Z, XU B, GAO B. An Image-segmentation-based Urban DTM Generation Method Using Airborne LiDAR Data[J]. IEEE Journal of Selected Topics in Applied Earth Observations & Remote Sensing, 2016, 9(1): 496-506.

[18] SITHOLE G, VOSSELMAN G. Filtering of Airborne Laser Scanner Data based on Segmented Point Clouds[J]. International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences, 2005, 36: 66-71.

[19] SHEN J,LIU J P,LIN X G,et al.Object-based Classification of Airborne Light Detection and Ranging Point Clouds in Human Settlements[J]. Sensor Letters, 2012, 10(1-2): 221-229.

[20] YAN Menglong,BLASCHKE T,LIU Yu,et al.An Object-based Analysis Filtering Algorithm for Airborne Laser Scanning[J]. International Journal of Remote Sensing, 2012, 33(22): 7099-7116.

[21] LIN Xiangguo, ZHANG Jixian. Segmentation-based Filtering of Airborne LiDAR Point Clouds By Progressive Densification of Terrain Segments[J]. Remote Sensing, 2014, 6(2): 1294-1326.

[22] XU S, VOSSELMAN G, ELBERINK O S. Multiple-Entity Based Classification of Airborne Laser Scanning Data in Urban Areas[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2014, 88: 1-15.

[23] RABBANI T. Automatic Reconstruction of Industrial Installations Using Point Clouds and Images[D]. Delft, the Netherlands:Netherlands Commission of Geodesy, 2006.

[24] ARYA S, MOUNT D M, NETANYAHU N S, et al. An Optimal Algorithm for Approximate Nearest Neighbor Searching Fixed Dimensions[J]. Journal of the ACM, 1998, 45(6): 891-923.

[25] 官云蘭, 程效軍, 施貴剛. 一種穩(wěn)健的點(diǎn)云數(shù)據(jù)平面擬合方法[J]. 同濟(jì)大學(xué)學(xué)報(自然科學(xué)版), 2008, 36(7): 981-984.

GUAN Yunlan, CHENG Xiaojun, SHI Guigang. A Robust Method for Fitting a Plane to Point Clouds[J]. Journal of Tongji University (Natural Science), 2008, 36(7): 981-984.

[26] EDELSBRUNNER H, KIRKPATRICK D, SEIDEL R. On the Shape of A Set of Points in the Plane[J]. IEEE Transactions on Information Theory, 1983, 29(4): 551-559.

[27] 亢曉琛, 劉紀(jì)平, 林祥國. 多核處理器的機(jī)載激光雷達(dá)點(diǎn)云并行三角網(wǎng)漸進(jìn)加密濾波方法[J]. 測繪學(xué)報, 2013, 42(3): 331-336. KANG Xiaochen, LIU Jiping, LIN Xiangguo. Parallel Filter of Progressive TIN Densification for Airborne LiDAR Point Cloud Using Multi-core CPU[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(3): 331-336.

(責(zé)任編輯：張艷玲)

Filtering of Point Clouds Using Fusion of Three Types of Primitives Including Points, Objects and Key Points

LIN Xiangguo，ZHANG Jixian，NING Xiaogang，DUAN Minyan，ZANG Yi

Chinese Academy of Surveying and Mapping, Beijing 100830, China

Primitive, being the basic processing unit, is one of the key factors to determine the accuracy and efficiency of point cloud filtering. Triangular irregular network (TIN) progressive densification (TPD) and object-based TIN progressive densification (OTPD) are two existing filtering methods, but single primitive is employed by them. A multiple-primitives-based TIN progressive densification (MPTPD) filtering method is proposed. It is composed of three key stages, including point cloud segmentation, extraction of key points of objects, the key-points-based judging of the objects. Specifically, point, object and the key points are the primitive of the above three stages respectively. Four testing datasets, including two airborne LiDAR and two photogrammetric point clouds, are used to verify the overall performances of the above three filtering methods. Experimental results suggest that the proposed MPTPD has the best overall performance. In the viewpoint of accuracy, MPTPD and OTPD have the similar accuracy. Moreover, compared with the TPD, MPTPD is able to reduce omission errors and total errors by 22.07% and 8.44% respectively. In the viewpoint of efficiency, under most of the cases, TPD is the highest, MPTPD is the second, and OTPD is the slowest. Moreover, the total time cost of MPTPD is only 57.93% of the one of OTPD.

filtering; LiDAR point cloud; photogrammetric point cloud; objects; triangular irregular network

The National Natural Science Foundations of China (No.41371405); The Foundation for Remote Sensing Young Talents by the National Remote Sensing Center of China;The Basic Research Fund of the Chinese Academy of Surveying and Mapping(No.777161103)

LIN Xiangguo(1981—), male, associate professor, post doctor, master supervisor, majors in image analysis and LiDAR data processing.

林祥國，張繼賢，寧曉剛，等.融合點(diǎn)、對象、關(guān)鍵點(diǎn)等3種基元的點(diǎn)云濾波方法[J].測繪學(xué)報，2016,45(11)：1308-1317.

10.11947/j.AGCS.2016.20160372.

LIN Xiangguo，ZHANG Jixian，NING Xiaogang，et al.Filtering of Point Clouds Using Fusion of Three Types of Primitives Including Points, Objects and Key Points[J]. Acta Geodaetica et Cartographica Sinica,2016,45(11):1308-1317. DOI:10.11947/j.AGCS.2016.20160372.

P237

A

1001-1595(2016)11-1308-10

國家自然科學(xué)基金(41371405)；遙感青年科技人才創(chuàng)新資助計劃；中國測繪科學(xué)研究院基本科研業(yè)務(wù)費(fèi)(777161103)

2016-07-29

修回日期： 2016-09-31

林祥國(1981—)，男，副研究員，博士后，碩士生導(dǎo)師，主要從事遙感影像分析、LiDAR數(shù)據(jù)處理方法研究。

E-mail： linxiangguo@casm.ac.cn