一種灰度體素結(jié)構(gòu)分割模型下的機(jī)載LiDAR3D濾波算法

　　摘要：針對(duì)現(xiàn)有的基于機(jī)載LiDAR數(shù)據(jù)的濾波算法未能充分利用數(shù)據(jù)提供的所有信息及其所采用的數(shù)據(jù)結(jié)構(gòu)表達(dá)復(fù)雜、存在信息損失等缺陷，提出了一種灰度體素結(jié)構(gòu)分割模型下的機(jī)載LiDAR3D濾波算法。算法首先以綜合利用機(jī)載LiDAR數(shù)據(jù)的高程及強(qiáng)度信息為目的將點(diǎn)云數(shù)據(jù)規(guī)則化為灰度(體素內(nèi)激光點(diǎn)的平均強(qiáng)度的離散化表示)體素結(jié)構(gòu)，然后基于各體素間的空間連通性和灰度相似性準(zhǔn)則，將灰度體素結(jié)構(gòu)分割并標(biāo)記為若干個(gè)3D連通區(qū)域，最后依據(jù)地面與其它目標(biāo)的高差特性提取與其對(duì)應(yīng)的3D連通區(qū)域。算法優(yōu)勢在于：基于體素結(jié)構(gòu)設(shè)計(jì)，為3D濾波算法;綜合利用了地面目標(biāo)的幾何及輻射特征，對(duì)比傳統(tǒng)濾波算法可應(yīng)用于更復(fù)雜的場景;濾波結(jié)果為3D地面體素形式，可直接用于創(chuàng)建地面3D模型。實(shí)驗(yàn)采用國際攝影測量與遙感協(xié)會(huì)(InternationalSocietyforPhotogrammetryandRemoteSensing,ISPRS)提供的不同密度的機(jī)載LiDAR基準(zhǔn)測試數(shù)據(jù)測試了鄰域尺度參數(shù)的敏感性及提出的算法的有效性，并和其他經(jīng)典濾波算法做對(duì)比。定量評(píng)價(jià)的結(jié)果表明，51鄰域?yàn)樽罴芽臻g鄰域尺度;點(diǎn)云密度為0.67點(diǎn)/m2的數(shù)據(jù)集1的濾波平均完整率、正確率及質(zhì)量分別為0.9611、0.9248及0.8934;點(diǎn)云密度為4點(diǎn)/m2的數(shù)據(jù)集2的濾波平均完整率、正確率及質(zhì)量分別為0.8490、0.8531及0.7404;對(duì)比其全經(jīng)典濾波算法本文算法在高密度點(diǎn)云數(shù)據(jù)濾波時(shí)表現(xiàn)更佳。

　　本文源自地球信息科學(xué)學(xué)報(bào),2020,22(11):2118-2127.《地球信息科學(xué)學(xué)報(bào)》(雙月刊)，1996年創(chuàng)刊，是綜合性學(xué)術(shù)刊物。《地球信息科學(xué)學(xué)報(bào)》融合了地理信息系統(tǒng)、全球定位系統(tǒng)、遙感、信息網(wǎng)絡(luò)等多學(xué)科的理論與技術(shù)，是綜合研究地球科學(xué)復(fù)雜系統(tǒng)的新領(lǐng)域。以地球系統(tǒng)信息流為主要研究對(duì)象;探討地球信息機(jī)理、地球信息認(rèn)知方法和地球信息時(shí)空圖譜技術(shù);實(shí)現(xiàn)對(duì)資源、環(huán)境協(xié)調(diào)發(fā)展，生態(tài)、區(qū)域持續(xù)發(fā)展等理論與實(shí)踐應(yīng)用的研究，應(yīng)對(duì)“數(shù)字地球”戰(zhàn)略與全球變化等科學(xué)問題。

　　1、引言

　　機(jī)載激光雷達(dá)(LightDetectionandRanging,LiDAR)數(shù)據(jù)濾波是目前點(diǎn)云數(shù)據(jù)處理領(lǐng)域中的關(guān)鍵技術(shù)，其精度直接影響數(shù)字高程模型(DigitalElevationModel,DEM)等產(chǎn)品的質(zhì)量。因此，濾波算法研究具有重要意義。

　　機(jī)載LiDAR系統(tǒng)可以快速獲取地表高精度、高密度的3D點(diǎn)云數(shù)據(jù)及其強(qiáng)度信息。但是，目前的濾波算法往往僅利用地面的幾何信息(如梯度[1,2]、連續(xù)[3,4,5,6]等)，這種數(shù)據(jù)利用方式未考慮點(diǎn)云強(qiáng)度信息在濾波中的輔助作用，未能充分發(fā)掘點(diǎn)云數(shù)據(jù)提供的所有信息，勢必影響濾波的準(zhǔn)確性。另外，已有濾波算法設(shè)計(jì)針對(duì)地形起伏不大、地形連續(xù)、地物較為單一的簡單場景，由此導(dǎo)致算法無法適用于地形起伏較大、地形不連續(xù)、地面附著物形態(tài)各異、規(guī)模不一、數(shù)量眾多的復(fù)雜場景[7]。

　　如果能綜合利用點(diǎn)云數(shù)據(jù)的幾何及強(qiáng)度信息并針對(duì)復(fù)雜地形特點(diǎn)設(shè)計(jì)算法，必將提高濾波結(jié)果的準(zhǔn)確性并促進(jìn)機(jī)載LiDAR在復(fù)雜場景下的濾波應(yīng)用。綜合利用所有信息的關(guān)鍵首先在于對(duì)機(jī)載LiDAR數(shù)據(jù)的表達(dá)方式。但傳統(tǒng)的數(shù)據(jù)結(jié)構(gòu)，如柵格格網(wǎng)[8,9]、不規(guī)則三角網(wǎng)(TriangulatedIrregularNetwork,TIN)[9]、點(diǎn)云[11]和八叉樹[12]應(yīng)用于機(jī)載LiDAR數(shù)據(jù)表達(dá)時(shí)均存在局限性：柵格格網(wǎng)和TIN為2.5D數(shù)據(jù)結(jié)構(gòu)，用其表達(dá)3D點(diǎn)云將導(dǎo)致信息損失并進(jìn)一步影響濾波結(jié)果的完整性;點(diǎn)云并未明晰表達(dá)激光點(diǎn)間的空間結(jié)構(gòu)及拓?fù)湫畔ⅲ纱藢?dǎo)致濾波算法設(shè)計(jì)困難、效率低下;八叉樹各節(jié)點(diǎn)的尺寸不一，各節(jié)點(diǎn)間的鄰接關(guān)系難以建立，這同樣增加了算法設(shè)計(jì)的難度。為了克服上述限制，采用一種更為簡單的真3D數(shù)據(jù)結(jié)構(gòu)表達(dá)機(jī)載LiDAR數(shù)據(jù)顯得尤為重要。同時(shí)，針對(duì)已有的濾波算法存在的無法適應(yīng)復(fù)雜場景的問題，有必要開發(fā)一個(gè)適應(yīng)范圍更廣泛的算法以更好地處理復(fù)雜多變場景。

　　因此，本文提出了一種基于灰度體素結(jié)構(gòu)分割模型下的機(jī)載LiDAR3D濾波算法。該算法首先將機(jī)載LiDAR數(shù)據(jù)規(guī)則化為灰度體素結(jié)構(gòu)，其中灰度為各體素內(nèi)激光點(diǎn)的平均強(qiáng)度值的離散化表示;然后基于體素間的空間連通性和灰度相似性準(zhǔn)則，將灰度體素結(jié)構(gòu)分割并標(biāo)記為若干個(gè)3D連通區(qū)域;最后依據(jù)地面的高程特性檢測地面對(duì)應(yīng)的3D連通區(qū)域。該算法的優(yōu)勢在于：(1)灰度體素結(jié)構(gòu)的真3D特性使其具備完備表達(dá)3DLiDAR點(diǎn)云的能力且該結(jié)構(gòu)同時(shí)融合了機(jī)載LiDAR數(shù)據(jù)的高程(各體素間的鄰接關(guān)系為高程接近的隱含表達(dá))和強(qiáng)度(各體素依據(jù)其內(nèi)的激光點(diǎn)強(qiáng)度賦值)信息，有助于綜合利用LiDAR數(shù)據(jù)的所有信息從而提高濾波精度;(2)3D濾波算法。利用地面目標(biāo)在灰度體素結(jié)構(gòu)中的“空間上具有連通性、灰度上表現(xiàn)出一致性”特性構(gòu)建濾波模型，基于灰度體素結(jié)構(gòu)分割模型、針對(duì)復(fù)雜場景中地面的幾何及輻射特征設(shè)計(jì)，明顯優(yōu)于僅利用地面幾何信息的濾波算法[1,2,3,4,5,6]的設(shè)計(jì)，可用于更復(fù)雜場景的濾波應(yīng)用。(3)濾波結(jié)果為3D地面體素形式，可直接用作3D地面模型。

　　2、基于體素分割的3D濾波

　　算法包括機(jī)載LiDAR數(shù)據(jù)規(guī)則化、灰度體素結(jié)構(gòu)分割及地面目標(biāo)檢測3個(gè)步驟。

　　2.1機(jī)載LiDAR數(shù)據(jù)規(guī)則化

　　規(guī)則化就是將包含場景目標(biāo)的空間范圍依據(jù)體素分辨率劃分成3D格網(wǎng)并依據(jù)格網(wǎng)單元(即體素)內(nèi)激光點(diǎn)的強(qiáng)度屬性對(duì)各個(gè)體素賦值為不同灰度等級(jí)的過程。其中，空間范圍由包圍機(jī)載LiDAR數(shù)據(jù)的軸向平行包圍盒確定;體素分辨率依據(jù)激光點(diǎn)間的平均點(diǎn)間距由式(1)確定;體素賦值方式如式(2)所示，將LiDAR數(shù)據(jù)中的各激光點(diǎn)映射到3D格網(wǎng)中，并對(duì)含有激光點(diǎn)的體素賦值激光點(diǎn)強(qiáng)度均值、不含有激光點(diǎn)的體素賦值0，進(jìn)而把上述體素值離散化到{0,1,???,255}得到體素灰度值。

　　式中：(Δx,Δy,Δz)為體素沿x、y和z方向的分辨率;Axy(Axz,Ayz)是LiDAR數(shù)據(jù)的XY(XZ,YZ)平面投影所得2D點(diǎn)集沿軸向最小外接矩形的面積;n為點(diǎn)數(shù)，式中取最小值是為了使構(gòu)建的體素結(jié)構(gòu)與LiDAR數(shù)據(jù)間的精度損失更少。

　　式中：i是激光點(diǎn)的索引;(xi,yi,zi)是第i個(gè)激光點(diǎn)的坐標(biāo);(ri,ci,li)是第i個(gè)體素的坐標(biāo)(行、列和層號(hào))。

　　另外，機(jī)載LiDAR數(shù)據(jù)中通常包含異常數(shù)據(jù)，如圖1所示。該數(shù)據(jù)源于低空飛行物、車輛等目標(biāo)對(duì)激光束的直接反射或樹干等對(duì)激光束的多次反射，為非真實(shí)目標(biāo)信息，其存在會(huì)影響體素結(jié)構(gòu)構(gòu)建的效率及準(zhǔn)確性，須在體素結(jié)構(gòu)構(gòu)建前予以剔除。剔除方案：統(tǒng)計(jì)LiDAR數(shù)據(jù)中各激光點(diǎn)高程的頻次，并通過直方圖可視化顯示;確定與真實(shí)目標(biāo)對(duì)應(yīng)的最高(低)高程閾值Th(Tl);對(duì)各激光點(diǎn)，若其高程高于Th或低于Tl，則判為異常點(diǎn)，予以剔除，否則保留。

　　圖1Samp31異常數(shù)據(jù)分布側(cè)視圖

　　2.2灰度體素結(jié)構(gòu)分割

　　分割的目標(biāo)就是將3D連通且灰度特征相似的體素合并到一個(gè)3D連通區(qū)域。假設(shè)2.1節(jié)所得灰度體素結(jié)構(gòu)中共包含l個(gè)3D連通區(qū)域，分割的任務(wù)就是將l個(gè)標(biāo)簽分配給體素結(jié)構(gòu)中的各個(gè)體素，以便屬于同一個(gè)3D連通區(qū)域的體素具有相同的標(biāo)簽，而屬于不同的3D連通區(qū)域的體素具有不同的標(biāo)簽。詳細(xì)方案：依次掃描灰度體素結(jié)構(gòu)中的非0值體素，直到掃描至第k個(gè)(k=1,2,…)未被標(biāo)記的體素。假設(shè)該體素的體素值為u、標(biāo)簽L1,…,Ld-1(其中，d是3D連通區(qū)域的索引，1≤d≤l)已被使用，選擇一個(gè)新的標(biāo)簽Ld。調(diào)用程序LABEL(k,u,Ld)(圖2)標(biāo)記與第k個(gè)體素3D連通(6、18、26或其它鄰域尺度的連通)且體素值相近(相似性準(zhǔn)則將在后文確定)的體素的標(biāo)簽為Ld。繼續(xù)掃描灰度體素結(jié)構(gòu)中的未被標(biāo)記的非0值體素，直到所有的體素都被標(biāo)記，得到若干個(gè)3D連通區(qū)域。

　　圖23D連通區(qū)域標(biāo)記LABEL(k,u,Ld)程序流程

　　在上述分割過程中應(yīng)用不同的鄰域尺度及體素值相似性準(zhǔn)則會(huì)得到不同的分割結(jié)果并由此影響后繼的濾波精度。為此，需面向應(yīng)用確定最優(yōu)值以保證地面的提取質(zhì)量。最佳鄰域尺度的確定見實(shí)驗(yàn)3.2節(jié)。最佳的體素值相似準(zhǔn)則方案如下(以實(shí)驗(yàn)數(shù)據(jù)Area3為例)：統(tǒng)計(jì)灰度體素結(jié)構(gòu)中的非0值體素的灰度值的頻率，并可視化顯示，如圖3所示。

　　圖3Area3的非0值灰度值頻率直方圖

　　由圖3可知：地物目標(biāo)的的灰度分布呈現(xiàn)多峰性(4個(gè)峰，分別位于1、28、80和133附近)且為多峰正態(tài)混合分布。則可借助高斯混合模型[13]對(duì)該分布進(jìn)行擬合，估計(jì)所得4個(gè)高斯分量的參數(shù)記做(μ2，σ2)，其中，j為高斯分布的索引，j=1、2、3、4。為了保證屬于同一正態(tài)分布的體素被分割到一個(gè)3D連通區(qū)域，可利用圖1中的峰值和谷值(位于10、61和103附近，分別記做v1,v2和v3，可利用局部極小值法確定)確定各個(gè)高斯分布的灰度范圍。以第二個(gè)高斯分布為例，若令μ2-m2l×σ2=v1，μ2+m2r×σ2=v2,則可確定m2l和m2r，由此可確定第二個(gè)高斯分布的灰度范圍[μ2-m2σ2，μ2+m2σ2]，其中，乘數(shù)m2=average{m2l,m2r}。同理可確定其它高斯分布的灰度范圍，并記做[rjl,rjr]。但是，相鄰的高斯分布間可能存在重疊。為了避免分布重疊且有利于后繼的濾波處理，各高斯分布的范圍分別設(shè)為(0,r2l),[r2l,r2r],(r2r,r3r],(r3r,r4r]。原因在于：第2個(gè)高斯分布對(duì)應(yīng)地面目標(biāo)(該先驗(yàn)知識(shí)可通過目視灰度體素模型獲取)，第1和第3個(gè)高斯分布的灰度范圍設(shè)置均以第2個(gè)為基準(zhǔn)設(shè)置以便準(zhǔn)確地獲取地面目標(biāo)。而第3和第4個(gè)高斯分布的分布重疊并未處理的原因在于算法目的在于分離地面目標(biāo)，此處即使植被和建筑物被混分也不影響算法的濾波精度。由此，前文所述的體素值相似性準(zhǔn)則確定轉(zhuǎn)換為了求取各個(gè)地物目標(biāo)的灰度范圍，圖2中的與u接近條件等價(jià)于“體素值位于第k個(gè)體素所對(duì)應(yīng)的目標(biāo)的統(tǒng)計(jì)灰度范圍內(nèi)”。

　　2.3地面目標(biāo)的3D連通區(qū)域檢測

　　利用地面低于周邊地物高度的特性，從分割所得的3D連通區(qū)域中分離地面目標(biāo)對(duì)應(yīng)的3D連通區(qū)域。詳細(xì)方案如下：

　　對(duì)任一3D連通區(qū)域，若其輪廓線的高程低于周圍地物一定的高度(高差閾值Te,Te=2m)，則作為候選，否則判定為非地面連通區(qū)域。其中，輪廓線的高程指的是位于輪廓線上的體素的平均高程。周圍地物的高程可用下述方案獲取：利用結(jié)構(gòu)元素[111;111;111]對(duì)該3D連通區(qū)域做膨脹，記膨脹處理所得外輪廓線的體素集合為Cs={va(ra,ca,la),a=1,2,…,q}，其中，s代表第s個(gè)3D連通區(qū)域，a是第s個(gè)3D連通區(qū)域的外輪廓線中各個(gè)體素的索引，q是第s個(gè)3D連通區(qū)域的外輪廓線包含的體素個(gè)數(shù)。va(ra,ca,la)∈Cs，搜尋與其平面坐標(biāo)相同的非0值體素，上述體素的平均高程即是第s個(gè)3D連通區(qū)域的周圍地物的平均高程。

　　計(jì)算各個(gè)候選3D連通區(qū)域的重心，若某一連通區(qū)域的重心高程值與其一定水平鄰域(如5×5)內(nèi)非0值體素的平均高程的高差小于等于某一閾值，則判定該3D連通區(qū)域?qū)?yīng)地面，否則對(duì)應(yīng)非地面。

　　3、實(shí)驗(yàn)數(shù)據(jù)、結(jié)果及討論

　　3.1實(shí)驗(yàn)數(shù)據(jù)

　　采用不同密度的兩組實(shí)測LiDAR點(diǎn)云數(shù)據(jù)作為實(shí)驗(yàn)數(shù)據(jù)檢驗(yàn)提出的算法的有效性和可行性。

　　數(shù)據(jù)集1：國際攝影測量與遙感協(xié)會(huì)(InternationalSocietyforPhotogrammetryandRemoteSensing,ISPRS)III/3工作組提供的濾波測試數(shù)據(jù)(https://www.itc.nl/isprs/wgIII-3/filtertest/downloadsites/)[14]:Samp11、12、21、22、23、24、31、41和42，圖4(a)以Samp31為例進(jìn)行展示。這些區(qū)域包括了不同場景中濾波可能遇到的主要困難，如粗差點(diǎn)的影響、復(fù)雜不規(guī)則形狀的大型建筑物、地物與地面相連、過街天橋、低矮植被等。點(diǎn)云數(shù)據(jù)記錄了首、末次回波及其強(qiáng)度信息，點(diǎn)云密度為0.67點(diǎn)/m2。另外，該數(shù)據(jù)已被手工分為地面和非地面點(diǎn)兩類，可用作標(biāo)準(zhǔn)數(shù)據(jù)定量評(píng)價(jià)提出的算法的精度。

　　圖4實(shí)驗(yàn)(點(diǎn)云)數(shù)據(jù)示例

　　數(shù)據(jù)集2:ISPRSIII/4工作組提供的目標(biāo)分類算法測試數(shù)據(jù)[15](數(shù)據(jù)需申請，http://www2.isprs.org/commissions/comm3/wg4/data-request-form2.html):Area2和Area3(網(wǎng)站未提供與Area1對(duì)應(yīng)的地面真實(shí)數(shù)據(jù)，無法進(jìn)行定量精度評(píng)價(jià)，因此該數(shù)據(jù)未被采用)，圖4(b)以Area3為例進(jìn)行展示。其中，Area2為被樹木圍繞的高層城市住宅;Area3為包含有獨(dú)立房屋和許多周圍樹木的住宅區(qū)。點(diǎn)云數(shù)據(jù)記錄了多次回波及其強(qiáng)度信息，點(diǎn)云密度為4點(diǎn)/m2。另外，該數(shù)據(jù)已被準(zhǔn)確的分類為瀝青地面、自然地面、建筑物、植被等類別，可用其中的地面真實(shí)數(shù)據(jù)作為參考數(shù)據(jù)定量評(píng)價(jià)提出的算法的精度。

　　3.2濾波過程及結(jié)果

　　實(shí)驗(yàn)以Samp31和Area3為例展示濾波過程及其結(jié)果，并對(duì)“鄰域尺度”參數(shù)進(jìn)行敏感性測試及分析。首先，Samp31和Area3經(jīng)異常數(shù)據(jù)剔除、空間范圍確定、體素分辨率計(jì)算、3D格網(wǎng)劃分、點(diǎn)云映射及體素賦值等步驟處理得到灰度體素結(jié)構(gòu)，如圖5所示。其中，各步驟的處理結(jié)果見表1。可目視圖5獲取各地物目標(biāo)的灰度分布情況，并作為先驗(yàn)知識(shí)用于隨后的體素值相似性準(zhǔn)則確定。

　　然后，對(duì)灰度體素結(jié)構(gòu)進(jìn)行分割。如2.2節(jié)所述，分割結(jié)果與各地物目標(biāo)的灰度范圍和鄰域尺度有關(guān)。其中，前者可由2.2節(jié)所述方案確定：由圖3可知第二個(gè)高斯分布對(duì)應(yīng)地面目標(biāo)，因此需首先確定其分布參數(shù)(μ2=33.34,σ2=12.79)，然后結(jié)合谷值10和61確定乘數(shù)m2l=1.82,m2r=2.16，并由此確定乘數(shù)m2=(1.82+2.16)/2=1.99，進(jìn)而確定第二個(gè)高斯分布的灰度范圍為[33.34-1.99×12.79,33.34+1.99×12.79]，即[8,59]。同理可確定其它高斯分布的灰度分布范圍。避免灰度分布重疊后的各個(gè)高斯分布的灰度范圍分別為(0,8),[8,59],(59,103],(103,255]。后者的最優(yōu)值實(shí)驗(yàn)確定方案如下：在相同條件下，不同的鄰域尺度分別被應(yīng)用到各個(gè)測試數(shù)據(jù)，對(duì)應(yīng)的濾波結(jié)果的誤差見表2。此處提出的算法所得地面數(shù)據(jù)是用體素表達(dá)的，而參考數(shù)據(jù)中則是離散的激光點(diǎn)。為和參考數(shù)據(jù)做對(duì)比進(jìn)而計(jì)算提出的算法的誤差，此處首先統(tǒng)計(jì)算法所得地面體素內(nèi)的激光點(diǎn)，然后和參考數(shù)據(jù)對(duì)比，進(jìn)而用I類誤差(將地面點(diǎn)錯(cuò)分為非地面點(diǎn)比例，Ie)、II類誤差(將非地面點(diǎn)錯(cuò)分為地面點(diǎn)比例，IIe)及總誤差(錯(cuò)分的地面點(diǎn)的比例，Te)等指標(biāo)定量描述算法誤差。

　　圖5構(gòu)建所得灰度體素模型

　　表1灰度體素模型構(gòu)建過程及結(jié)果

　　由表2可知：(1)數(shù)據(jù)集1中，采用6、18、26、51及56鄰域時(shí)，提出的算法的平均總誤差分別為0.0907、0.0765、0.0723、0.0685及0.0800。也即，從總誤差指標(biāo)來看，點(diǎn)云密度為0.67點(diǎn)/m2時(shí)，51鄰域是最佳鄰域尺度。(2)數(shù)據(jù)集2中，采用6、18、26、51及56鄰域時(shí)，提出的算法的平均總誤差分別為0.1029、0.1025、0.1019、0.0838及0.0962。也即，從總誤差指標(biāo)來看，點(diǎn)云密度為4點(diǎn)/m2時(shí)，51鄰域是最佳鄰域尺度。(3)鄰域尺度的增加并不意味著總誤差的必然降低。提出的算法的思想是：地面目標(biāo)信息可以通過在灰度體素模型中定義的連通性和灰度相似度來傳播。以6鄰域?yàn)槔孛婺繕?biāo)信息的傳播只能由中心體素向上、下或4個(gè)基本方向傳播。由此導(dǎo)致采用6鄰域只能將平坦地區(qū)的地面體素合并到一個(gè)3D連通區(qū)域并繼而被正確分離。如果能增加鄰域尺度，如18、26鄰域，其傳播方向的增加可能會(huì)納入更多的地面體素，從而降低濾波的I類誤差(如表2中采用18、26鄰域時(shí)的Ie明顯低于采用6鄰域算法的Ie)。但是，如果鄰域尺度過大，則可能將一些非地面體素錯(cuò)分為地面體素，從而導(dǎo)致II類誤差增加從而降低濾波的準(zhǔn)確性(如表2中56鄰域算法的IIe明顯高于51鄰域算法的IIe)。

　　51鄰域尺度下的體素模型分割結(jié)果見圖6。由圖7可知：地面目標(biāo)被分割為多個(gè)獨(dú)立的3D連通區(qū)域。

　　最后，基于各個(gè)3D連通區(qū)域的特性分離地面目標(biāo)。Samp11、12、21、22、23、24、31、41、42、Area2和Area3所得地面體素?cái)?shù)分別為18517、18728、7095、15854、11253、3893、12366、3096、9905、39729、39154,Samp31和Area3所得地面體素見圖7。

　　表2提出的算法在不同鄰域尺度下的的誤差統(tǒng)計(jì)

　　圖6灰度體素模型分割結(jié)果(各連通區(qū)域用不同顏色表示)

　　圖7提出的濾波算法的濾波結(jié)果

　　3.3濾波算法定量精度評(píng)價(jià)

　　為了定量評(píng)價(jià)提出的算法精度，以基于像元的評(píng)價(jià)方法[16]進(jìn)一步統(tǒng)計(jì)了濾波結(jié)果的完整率(Rcom)、正確率(Rcor)、質(zhì)量(Rq)及Kappa系數(shù)，如表3所示。

　　由表3可知，數(shù)據(jù)集1的平均完整率、正確率、質(zhì)量及Kappa系數(shù)分別為0.9611、0.9248、0.8934及0.8484;數(shù)據(jù)集2的平均完整率、正確率、質(zhì)量及Kappa系數(shù)分別為0.8490、0.8531、0.7404及0.7922。為了探究影響算法完整率和正確率的因素，考察了濾波結(jié)果的誤差分布情況，如圖8所示。

　　結(jié)合圖8分析可知，影響提出的算法正確率的主要因素包括：某些和地面目標(biāo)相鄰且強(qiáng)度接近的非地面目標(biāo)被錯(cuò)分為地面目標(biāo);算法中與地面對(duì)應(yīng)的連通區(qū)域分離僅利用了地面低于周圍地物的特性，可能導(dǎo)致錯(cuò)分。影響提出的算法完整率的因素包括：某些高程突變的地面被錯(cuò)分。

　　表3基于像元的評(píng)價(jià)方式下的提出算法的算法精度

　　圖8誤差分布頂視圖

　　3.4濾波算法精度對(duì)比

　　針對(duì)數(shù)據(jù)集1，提出的算法同已有經(jīng)典濾波算法的總誤差對(duì)比見表4。其中，Axelsson[10]為2003年前發(fā)表的精度最高的濾波算法，其他則為2003年后發(fā)表的在多個(gè)測試數(shù)據(jù)中誤差低于Axelsson[10]的算法。

　　由表4可知，提出的算法對(duì)比Axelsson[10]改進(jìn)了9個(gè)樣本中的3個(gè)樣本精度;提出的算法對(duì)比已有的精度最高的3D濾波算法[3]精度相當(dāng)，改進(jìn)了9個(gè)樣本中的2個(gè)樣本精度。

　　鑒于目前利用數(shù)據(jù)集二的均為分類算法(分為建筑物、樹、低矮植被、天然地面(naturalground)及瀝青地面(asphaltground)等類)，為了將提出的算法同上述算法(參見鏈接，Area2和Area3的鏈接分別為http://www2.isprs.org/commissions/comm3/wg4/results/a2_detect.html和http://www2.isprs.org/commissions/comm3/wg4/results/a3_detect.html)做對(duì)比，本文在提出的算法基礎(chǔ)上進(jìn)一步利用天然地面的激光反射強(qiáng)度值大于瀝青地面的激光反射強(qiáng)度值這一特性對(duì)二者加以區(qū)分，進(jìn)而統(tǒng)計(jì)上述檢測的完整率、正確率及質(zhì)量并和已有算法做對(duì)比，如表5所示。其中，天然和瀝青地面的區(qū)分細(xì)則為：若某一已判做地面的3D連通區(qū)域的灰度(3D連通區(qū)域內(nèi)各個(gè)體素的強(qiáng)度均值)大于某閾值，則為天然地面，否則為瀝青地面。閾值的選取規(guī)則為：統(tǒng)計(jì)已判做地面的各3D連通區(qū)域的灰度頻率，設(shè)置頻率直方圖中的谷值作為灰度閾值Ti。

　　由表5可知，提出的算法對(duì)比其他已有分類算法，其質(zhì)量指標(biāo)均為最高值，這也說明了算法的有效性。

　　綜合上述算法精度對(duì)比實(shí)驗(yàn)可知，提出的算法在低密度點(diǎn)云數(shù)據(jù)濾波中與已有的3D濾波算法精度相當(dāng)，但隨著點(diǎn)云密度的增加，其濾波精度明顯高于其它濾波算法，優(yōu)勢明顯。

　　表4提出的算法和其他算法的總誤差對(duì)比

　　表5提出的算法與其他算法的自然地面分類精度對(duì)比

　　3.5討論

　　提出的算法的濾波結(jié)果由輸入?yún)?shù)如閾值(Te、Ti)、目標(biāo)的統(tǒng)計(jì)灰度范圍和鄰域尺度等參數(shù)決定。其中，Ti和目標(biāo)的統(tǒng)計(jì)灰度范圍是根據(jù)實(shí)際數(shù)據(jù)源設(shè)置的，使用本文給出的設(shè)置方案，可以很容易地確定“數(shù)據(jù)源”型閾值，并不能限制提出的算法的普適性。Te的設(shè)定是經(jīng)驗(yàn)型的，因?yàn)榈孛嫱ǔ１绕渲車慕ㄖ镏辽俚?m，考慮到地面周圍還有低矮植被，所以本文設(shè)置高度差閾值為2m。點(diǎn)云密度約為0.67和4點(diǎn)/m2時(shí)的鄰域尺度可以直接使用51鄰域，因?yàn)楸疚囊炎C明51鄰域?yàn)闉V波的最佳鄰域尺度。因此，本文提出的算法與地面的形態(tài)、坡度及連續(xù)性無關(guān)，可適用于不同復(fù)雜度場景的地面目標(biāo)提取，具有普遍性。

　　4、結(jié)論

　　本文面向機(jī)載LiDAR數(shù)據(jù)提出了一種3D濾波算法。算法首先將機(jī)載LiDAR數(shù)據(jù)規(guī)則化為灰度體素結(jié)構(gòu)，然后利用體素間的連通性和灰度相似性將其分割成多個(gè)3D連通區(qū)域，最后根據(jù)地面和其他目標(biāo)間存在高差這一特性提取與地面對(duì)應(yīng)的3D連通區(qū)域。該算法基于灰度體素結(jié)構(gòu)設(shè)計(jì)，以3D連通區(qū)域分割理論為基礎(chǔ)，綜合利用了LiDAR點(diǎn)云數(shù)據(jù)包含的高程及強(qiáng)度信息。對(duì)比僅利用幾何信息的濾波算法，其將強(qiáng)度信息用于濾波輔助決策，從而為地面和非地面目標(biāo)的區(qū)分提供了更有效的信息，因而可適用于更復(fù)雜的場景。實(shí)驗(yàn)采用不同密度的ISPRS基準(zhǔn)測試數(shù)據(jù)測試了鄰域尺度參數(shù)的敏感性及提出的濾波算法的有效性，并和其他經(jīng)典的濾波算法的精度進(jìn)行對(duì)比。定量評(píng)價(jià)的結(jié)果表明，51鄰域?yàn)樽罴芽臻g鄰域尺度;點(diǎn)云密度為0.67(4)點(diǎn)/m2的數(shù)據(jù)集1(2)的濾波的平均完整率、正確率及質(zhì)量分別為0.9611(0.8490)、0.9248(0.8531)及0.8934(0.7404);提出的算法相對(duì)其他經(jīng)典濾波算法在高密度點(diǎn)云數(shù)據(jù)濾波時(shí)表現(xiàn)更佳。總體而言，提出的算法有助于綜合利用幾何及輻射信息以提高復(fù)雜場景下地面提取結(jié)果的準(zhǔn)確性。但是，算法僅支持高程或強(qiáng)度不同的地面目標(biāo)的有效分離。為支持更復(fù)雜場景的地面目標(biāo)分類，后續(xù)研究可嘗試構(gòu)建點(diǎn)云數(shù)據(jù)的多值體素結(jié)構(gòu)模型(即體素賦值體素內(nèi)激光點(diǎn)的色彩信息，可通過LiDAR和光學(xué)影像融合數(shù)據(jù)或者多光譜LiDAR數(shù)據(jù)等方式獲取激光點(diǎn)的色彩信息)并在此基礎(chǔ)上提出相應(yīng)的地面目標(biāo)分離算法。

　　參考文獻(xiàn)：

　　[5]劉凱斯,王彥兵,宮輝力,等.機(jī)載LiDAR點(diǎn)云數(shù)據(jù)的二面角濾波算法[J].地球信息科學(xué)學(xué)報(bào),2018,20(4):414-421.

　　[6]王麗英,王圣,徐艷,等.結(jié)合體元數(shù)據(jù)結(jié)構(gòu)的機(jī)載LIDAR建筑物檢測[J].中國圖象圖形學(xué)報(bào),2017,22(10):1436-1446.

　　[7]張杰.機(jī)載LiDAR點(diǎn)云濾波及特征提取技術(shù)研究[D].徐州:中國礦業(yè)大學(xué),2017.

　　[8]孫美玲.機(jī)載LiDAR數(shù)據(jù)濾波及城區(qū)汽車目標(biāo)檢測方法研究[D].成都:西南交通大學(xué),2014.

　　[9]孫蒙,顧和和.基于微分形態(tài)學(xué)斷面的機(jī)載LiDAR數(shù)據(jù)濾波新方法[J].大地測量與地球動(dòng)力學(xué),2016,36(7):591-594,599.

　　[12]周曉明.機(jī)載激光雷達(dá)點(diǎn)云數(shù)據(jù)濾波算法的研究與應(yīng)用[D].鄭州:解放軍信息工程大學(xué),2011.

　　[13]趙泉華,石雪,王玉,等.可變類空間約束高斯混合模型遙感圖像分割[J].通信學(xué)報(bào),2017,38(2):34-43.