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

林祥國(guó)，張繼賢，寧曉剛，段敏燕，臧 藝

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

?

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

林祥國(guó)，張繼賢，寧曉剛，段敏燕，臧 藝

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

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

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

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

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

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

1 本文的濾波方法

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

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

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

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

第1步，估計(jì)法向量和殘差，處理過(guò)程如下。

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

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

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

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

圖2 表面生長(zhǎng)過(guò)程的示意圖Fig.2 Illustration of the surface growing process

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

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

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

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

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

(6) 結(jié)束。

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

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

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

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

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

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

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

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

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

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

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

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

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

1.3 基于關(guān)鍵點(diǎn)的對(duì)象類(lèi)別判別

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

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

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

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

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

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

第5步：迭代的判別對(duì)象類(lèi)別。子步驟包括：

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

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

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

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

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

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

2 試驗(yàn)與分析

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

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

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

對(duì)4個(gè)試驗(yàn)數(shù)據(jù)進(jìn)行濾波時(shí)，使用的相關(guān)參數(shù)的值見(jiàn)表2。其中，TPD、OTPD和MPTPD3種方法對(duì)1.3節(jié)的參數(shù)b、θ、d、t采用了相同的參數(shù)值，OTPD和MPTPD兩種方法對(duì)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 精度評(píng)價(jià)

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

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

表3中，TPD方法的Ⅰ、Ⅱ、T3類(lèi)誤差的平均值分別為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ù)字說(shuō)明OTPD方法與MPTPD方法的各類(lèi)誤差均十分接近。且，一類(lèi)誤差I(lǐng)和總誤差T均呈現(xiàn)TPD>OTPD≈MPTPD的趨勢(shì)，而二類(lèi)誤差I(lǐng)I呈現(xiàn)TPD

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

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

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

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

試驗(yàn)數(shù)據(jù)誤差類(lèi)型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 效率評(píng)價(jià)

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

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

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

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

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

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

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

2.4 分析與討論

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

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

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

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

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

3 結(jié) 論

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

[1] 黃先鋒, 李卉, 王瀟, 等. 機(jī)載LiDAR數(shù)據(jù)濾波方法評(píng)述[J]. 測(cè)繪學(xué)報(bào), 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ìng)雪, 朱慶, 王偉璽. 多匹配基元集成的多視影像密集匹配方法[J]. 測(cè)繪學(xué)報(bào), 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á)測(cè)量技術(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é)報(bào)(信息科學(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é)報(bào), 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]. 測(cè)繪學(xué)報(bào), 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é)報(bào)(自然科學(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ì)平, 林祥國(guó). 多核處理器的機(jī)載激光雷達(dá)點(diǎn)云并行三角網(wǎng)漸進(jìn)加密濾波方法[J]. 測(cè)繪學(xué)報(bào), 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.

林祥國(guó)，張繼賢，寧曉剛，等.融合點(diǎn)、對(duì)象、關(guān)鍵點(diǎn)等3種基元的點(diǎn)云濾波方法[J].測(cè)繪學(xué)報(bào)，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

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

2016-07-29

修回日期： 2016-09-31

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

E-mail： linxiangguo@casm.ac.cn