一種灰度體素結構分割模型下的機載LiDAR3D濾波算法

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

　　本文源自地球信息科學學報,2020,22(11):2118-2127.《地球信息科學學報》(雙月刊)，1996年創刊，是綜合性學術刊物。《地球信息科學學報》融合了地理信息系統、全球定位系統、遙感、信息網絡等多學科的理論與技術，是綜合研究地球科學復雜系統的新領域。以地球系統信息流為主要研究對象;探討地球信息機理、地球信息認知方法和地球信息時空圖譜技術;實現對資源、環境協調發展，生態、區域持續發展等理論與實踐應用的研究，應對“數字地球”戰略與全球變化等科學問題。

　　1、引言

　　機載激光雷達(LightDetectionandRanging,LiDAR)數據濾波是目前點云數據處理領域中的關鍵技術，其精度直接影響數字高程模型(DigitalElevationModel,DEM)等產品的質量。因此，濾波算法研究具有重要意義。

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

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

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

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

　　算法包括機載LiDAR數據規則化、灰度體素結構分割及地面目標檢測3個步驟。

　　2.1機載LiDAR數據規則化

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

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

　　式中：i是激光點的索引;(xi,yi,zi)是第i個激光點的坐標;(ri,ci,li)是第i個體素的坐標(行、列和層號)。

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

　　圖1Samp31異常數據分布側視圖

　　2.2灰度體素結構分割

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

　　圖23D連通區域標記LABEL(k,u,Ld)程序流程

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

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

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

　　2.3地面目標的3D連通區域檢測

　　利用地面低于周邊地物高度的特性，從分割所得的3D連通區域中分離地面目標對應的3D連通區域。詳細方案如下：

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

　　計算各個候選3D連通區域的重心，若某一連通區域的重心高程值與其一定水平鄰域(如5×5)內非0值體素的平均高程的高差小于等于某一閾值，則判定該3D連通區域對應地面，否則對應非地面。

　　3、實驗數據、結果及討論

　　3.1實驗數據

　　采用不同密度的兩組實測LiDAR點云數據作為實驗數據檢驗提出的算法的有效性和可行性。

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

　　圖4實驗(點云)數據示例

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

　　3.2濾波過程及結果

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

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

　　圖5構建所得灰度體素模型

　　表1灰度體素模型構建過程及結果

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

　　51鄰域尺度下的體素模型分割結果見圖6。由圖7可知：地面目標被分割為多個獨立的3D連通區域。

　　最后，基于各個3D連通區域的特性分離地面目標。Samp11、12、21、22、23、24、31、41、42、Area2和Area3所得地面體素數分別為18517、18728、7095、15854、11253、3893、12366、3096、9905、39729、39154,Samp31和Area3所得地面體素見圖7。

　　表2提出的算法在不同鄰域尺度下的的誤差統計

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

　　圖7提出的濾波算法的濾波結果

　　3.3濾波算法定量精度評價

　　為了定量評價提出的算法精度，以基于像元的評價方法[16]進一步統計了濾波結果的完整率(Rcom)、正確率(Rcor)、質量(Rq)及Kappa系數，如表3所示。

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

　　結合圖8分析可知，影響提出的算法正確率的主要因素包括：某些和地面目標相鄰且強度接近的非地面目標被錯分為地面目標;算法中與地面對應的連通區域分離僅利用了地面低于周圍地物的特性，可能導致錯分。影響提出的算法完整率的因素包括：某些高程突變的地面被錯分。

　　表3基于像元的評價方式下的提出算法的算法精度

　　圖8誤差分布頂視圖

　　3.4濾波算法精度對比

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

　　由表4可知，提出的算法對比Axelsson[10]改進了9個樣本中的3個樣本精度;提出的算法對比已有的精度最高的3D濾波算法[3]精度相當，改進了9個樣本中的2個樣本精度。

　　鑒于目前利用數據集二的均為分類算法(分為建筑物、樹、低矮植被、天然地面(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)做對比，本文在提出的算法基礎上進一步利用天然地面的激光反射強度值大于瀝青地面的激光反射強度值這一特性對二者加以區分，進而統計上述檢測的完整率、正確率及質量并和已有算法做對比，如表5所示。其中，天然和瀝青地面的區分細則為：若某一已判做地面的3D連通區域的灰度(3D連通區域內各個體素的強度均值)大于某閾值，則為天然地面，否則為瀝青地面。閾值的選取規則為：統計已判做地面的各3D連通區域的灰度頻率，設置頻率直方圖中的谷值作為灰度閾值Ti。

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

　　綜合上述算法精度對比實驗可知，提出的算法在低密度點云數據濾波中與已有的3D濾波算法精度相當，但隨著點云密度的增加，其濾波精度明顯高于其它濾波算法，優勢明顯。

　　表4提出的算法和其他算法的總誤差對比

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

　　3.5討論

　　提出的算法的濾波結果由輸入參數如閾值(Te、Ti)、目標的統計灰度范圍和鄰域尺度等參數決定。其中，Ti和目標的統計灰度范圍是根據實際數據源設置的，使用本文給出的設置方案，可以很容易地確定“數據源”型閾值，并不能限制提出的算法的普適性。Te的設定是經驗型的，因為地面通常比其周圍的建筑物至少低3m，考慮到地面周圍還有低矮植被，所以本文設置高度差閾值為2m。點云密度約為0.67和4點/m2時的鄰域尺度可以直接使用51鄰域，因為本文已證明51鄰域為濾波的最佳鄰域尺度。因此，本文提出的算法與地面的形態、坡度及連續性無關，可適用于不同復雜度場景的地面目標提取，具有普遍性。

　　4、結論

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

　　參考文獻：

　　[5]劉凱斯,王彥兵,宮輝力,等.機載LiDAR點云數據的二面角濾波算法[J].地球信息科學學報,2018,20(4):414-421.

　　[6]王麗英,王圣,徐艷,等.結合體元數據結構的機載LIDAR建筑物檢測[J].中國圖象圖形學報,2017,22(10):1436-1446.

　　[7]張杰.機載LiDAR點云濾波及特征提取技術研究[D].徐州:中國礦業大學,2017.

　　[8]孫美玲.機載LiDAR數據濾波及城區汽車目標檢測方法研究[D].成都:西南交通大學,2014.

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

　　[12]周曉明.機載激光雷達點云數據濾波算法的研究與應用[D].鄭州:解放軍信息工程大學,2011.

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