陳曦亮,王 雪,畢曉偉

(1.長安大學地質工程與測繪學院，陜西 西安 710064; 2.咸陽師范學院,陜西 咸陽 712000)

0 引 言

LiDAR(Light Detection and Ranging)技術是一種主動式的遙感技術，這一技術使用的是激光設備發(fā)出的脈沖信號照射到地面，記錄每個脈沖信號從發(fā)出到打到地面點上返回所用的時間，利用時間和光速就可以計算出傳感器設備到地面某點的距離，運用GPS全球定位系統(tǒng)和INS慣性導航系統(tǒng)確認激光發(fā)射中心的三維坐標和姿態(tài)傾角[1]。獲得的點云數(shù)據(jù)有目標地面點的三維坐標信息，可以依次通過后續(xù)處理得到相關產(chǎn)品。這項技術具有可以全天候作業(yè)，可以獲取常規(guī)攝影測量無法獲取的困難地區(qū)的三維信息，作業(yè)周期短，生產(chǎn)費用相對較低，不需要昂貴的測圖設備，在困難地區(qū)精度高于攝影測量的精度以及可以直接獲取目標地物的DEM而不需要再次進行新的人工處理等優(yōu)點[2]。

在對獲取的點云數(shù)據(jù)進行處理的過程中，濾波的過程是相當關鍵的一步，如何找到一個最好的濾波方法是很多人研究的目的[3-4]。濾波方法有很多，具體可分為數(shù)學形態(tài)學濾波、TIN濾波、統(tǒng)計分析濾波、基于調整點云表面的濾波、分類與分割的濾波以及多尺度濾波等[5]，還有布料模擬濾波。每種方法都有其優(yōu)缺點?；诜诸惡头指畹姆椒ū容^適合城市地區(qū)，但對植被茂盛的地區(qū)有效性較低。TIN濾波的可靠性得到過很多證實，但對于地形不連續(xù)的區(qū)域濾波效果不太好[6]。統(tǒng)計分析濾波可以避免人為的參數(shù)調節(jié)，但穩(wěn)定性較差。多尺度濾波借助多種策略對非地面點進行識別整合，對于平臺地區(qū)效果較好，但對地形較復雜地區(qū)就不太適合。

2003年Zhang等[7]提出了一種漸進式的數(shù)學形態(tài)學濾波，該算法通過逐步迭代增大濾波窗口，并且根據(jù)坡度值動態(tài)變化高差閾值，但其地形自適應性較弱，設置高差閾值的時候也不能很好地顧及地形特點。2015年隋立春等[8]提出了一種改進的數(shù)學形態(tài)學濾波算法，提出了在“開”算子基礎上增加一個帶寬參數(shù)，提高了地形自適應性。2016年Zhang等[9]提出了通過計算機布料模擬的濾波算法。2018年，張昌賽等[10]對布料模擬的算法進行了適用性分析，給出了總結：布料模擬濾波算法具有效率高，在地形平坦地區(qū)精度也高，但在一些地形復雜的區(qū)域效果就欠缺一些。2020年，鄭津津等[11]基于布料模擬濾波和隨機森林進行了點云分類過程，并且提出了基于地面點進行高程歸一化的思路。本文圍繞這2種濾波算法來進行改進。利用這2類算法的特性來改進算法，使地面點和非地面點可以更高精度地分離，并適應多種地形。

1 數(shù)據(jù)預處理

在對點云數(shù)據(jù)進行濾波之前，需要對點云數(shù)據(jù)先進行一些預處理步驟，包括空間索引建立和一些粗差噪聲點的剔除。

1.1 空間索引的建立

由于點云數(shù)據(jù)本身十分龐大，因此在數(shù)據(jù)處理前將數(shù)據(jù)本身進行一定的組織化是十分必要的步驟，本文采用空間格網(wǎng)索引結構來組織點云數(shù)據(jù)，本方法在2004年由Cho等[12]將其引入到 LiDAR 點云數(shù)據(jù)預先處理中。其基本原理是將離散點云所分布的區(qū)域用一個個大小相同的虛擬格網(wǎng)分開，每個格網(wǎng)就相當于點云的一個子空間容器，記錄落入其中的點云子集。對于激光點云創(chuàng)建的虛擬格網(wǎng)如圖1所示，其中，左圖顯示的是立體結構，黑點表示離散的三維激光點云各個點，灰色的立方體則表示所建立的其中之一的虛擬格網(wǎng)；右圖則顯示的是平面結構。白色小圓圈表示離散的激光腳點，正方形表示創(chuàng)建的虛擬格網(wǎng)。

圖1 虛擬格網(wǎng)的空間索引

1.2 噪聲點的剔除

點云數(shù)據(jù)中常見的粗差(噪聲點)包括高位粗差和低位粗差，如圖2所示，機載的LiDAR系統(tǒng)在飛行中受到了低空飛行物(如飛機、鳥類等)的影響，誤將這些物體反射回來的信號當作了被測物體，高位粗差同時又叫做空中噪聲點。低位粗差則是由于多路徑因素影響又或者是激光測距儀自身的誤差帶來的極低點誤差，低位粗差一般是低于實際地形的噪聲粗差點，又稱低點[13]。

去噪算法的原理：該算法對每個點的領域進行搜索，計算從一個區(qū)域點到周圍相鄰的區(qū)域點的距離，計算該區(qū)域點的平均距離和標準差，如果一個區(qū)域的距離平均值大于最大距離(最大距離=中值+標準差×標準差倍數(shù))，則將該點視為噪聲的點，通過選擇適當?shù)膮?shù)，如果某點和周圍領域點的平均距離大于最大距離，可以將其視為粗差噪聲點。消除在點云中的粗差[14]。

2 算法原理

2.1 漸進式數(shù)學形態(tài)學濾波

數(shù)學形態(tài)學理論用于LiDAR濾波的時候，有2種算子，即“開”算子和“閉”算子?！伴_”算子先進行腐蝕處理，再進行膨脹處理，“閉”算子則正好相反。本文使用開算子。腐蝕和膨脹操作如圖3和圖4所示。結構B為窗口模板[15]。

圖3 腐蝕操作

圖4 膨脹操作

“開”算子的算法流程如下：

1)腐蝕處理。掃描搜索全部的LiDAR點云數(shù)據(jù),以其中的一點為中心開一個W×W大小的窗口(搜索模板)，該窗口在整個的點云數(shù)據(jù)圖像上進行逐行逐列的移動，這時記錄各點的高程值并且進行比較，找到里面最小的那個值，使用這個最小的高程值代替掃描窗口中心點的高程值，直至所有點的高程都腐蝕完后才結束腐蝕步驟。

2)膨脹處理。再次掃描搜索所有的點云數(shù)據(jù),使用和腐蝕操作同樣大小的搜索窗口,并使該窗口進行逐行逐列的移動。此時,各個點的高程值均已被腐蝕操作時窗口內最小的高程值所代替,同樣比較窗口各點高程值的大小，用窗口內高程值最大的值代替窗口中心點的高程值。繼續(xù)循環(huán)該操作，直到所有點的高程值更新完畢。

3)地面點判定和提取。假設P點操作之前的高程值為Zp，dhT為要設置的閾值。開算子的操作結束后，對每一個激光腳點是否為地面點進行判斷。其依據(jù)是：如果P點膨脹操作之后的高程值和其原始高程值Zp之差的絕對值小于或等于閾值dhT，則判定P點為地面點，若大于閾值就判定為非地面點[16]。

這里采用漸進式改變窗口大小的方法，其窗口大小隨著迭代步長的增長而指數(shù)增長，關系式為Wk=2×pow(b,k)+1，其中，Wk為現(xiàn)窗口大小，b為設置的初始窗口大小，k為迭代次數(shù)。窗口大小增長到設定的最大窗口值為止。

利用研究區(qū)域的坡度值來進行高差閾值的計算。高差閾值dhT和坡度值s以及窗口大小Wk的關系如式(1)所示：

(1)

其中，C值表示單元格網(wǎng)的大小值。

2.2 布料模擬濾波

布料模擬濾波(Cloth Simulation Filtering, CSF)算法的過程可以看作是一種簡單的模擬布料的物理過程。該算法假設在地表區(qū)域的正上方存在著一塊虛擬布料，在重力的作用下該布料最終會覆蓋在下方的地形表面上。假設這塊虛擬布料足夠柔軟，則布料落下來的最終形狀剛好便是對應下方區(qū)域點云的數(shù)字表面模型(DSM)。同理，如果把原始地形反向翻轉過來，同時布料具有一定的硬度,則同樣在重力作用下布料的最終形狀就是點云的數(shù)字高程模型(DEM)[17]。布料模擬點云濾波算法的原理如圖5所示。

圖5 布料模擬濾波示意圖

計算機里面所模擬的布料是由大量的離散節(jié)點構成的一個格網(wǎng)，格網(wǎng)中的節(jié)點相互聯(lián)系并由連接線進行連接，同時假設各格網(wǎng)點均為具有恒定質量且無尺寸的質點，因此這種格網(wǎng)的方式也稱作為質點彈簧模型[18]。要模擬布料在自由落體中某一時刻的形狀，需要計算出各布料點的空間位置。假設各布料節(jié)點只在垂直方向進行移動，而且只受重力與鄰近節(jié)點的相互作用力共2種力。當只考慮重力作用時，由牛頓第二定律可知，布料點的空間位置和其所受的作用力之間的關系可由式(2)確定：

X(t+△t)=2X(t)-X(t-△t)+G△t2/m

(2)

其中，m表示布料點的質量，G表示重力加速度常數(shù)，X(t)表示某一時刻的格網(wǎng)節(jié)點位置，△t表示時間步長。當時間步長和格網(wǎng)節(jié)點初始位置已知時，可以解算出格網(wǎng)節(jié)點的當前位置。

為了約束粒子在反轉表面空白區(qū)域的反轉問題，還需考慮粒子內部因素導致的位移 。任意選取2個相鄰的粒子，如果2個粒子都是可移動的，則令二者往相反的方向移動同樣的距離；如果一個是不可移動的，則移動另一個；如果兩者具有相同的高度，則不進行移動[19]。位移量可以通過式(3)進行計算：

(3)

其中，d為粒子的位移量；當粒子可移動時，b等于1，不可移動時b為0；pi為p0的相鄰粒子，n是點進行標準化到垂直方向上的單位向量(0，0，1)T。

布料模擬濾波算法的詳細步驟如下：

1)上下翻轉原始點云，對點云空間坐標進行轉換。

2)設定格網(wǎng)分辨率，將原始點云進行虛擬格網(wǎng)化。把所有激光點與格網(wǎng)上的粒子在同一個水平面投影，找到每一個粒子的最鄰近點，記錄其投影前的高程。

3)對于每一個可移動的格網(wǎng)粒子，計算其受到重力影響產(chǎn)生的位移，并將該位置與其對應激光腳點的高程進行比較，若布料格網(wǎng)節(jié)點高程小于或等于激光腳點高程，則將該節(jié)點移動到激光腳點的位置，并標記為不可移動點。

4)計算各布料格網(wǎng)節(jié)點受鄰近點作用力后需要移動的位移量。

5)重復步驟3和步驟4，當各節(jié)點的高程變化開始變得足夠小的時候或者迭代次數(shù)達到最大，布料模擬過程就可以結束。

6)區(qū)分地面點與非地面點，如果LiDAR點與模擬粒子之間的距離小于預先設置的閾值，則認為其是地面點，反之則認為其為非地面點。

2.3 算法的改進

2.3.1 2種傳統(tǒng)方法存在的問題

數(shù)學形態(tài)學算法的難點在于窗口大小的選取，窗口大小的參數(shù)決定了濾波的細節(jié)程度。窗口太大時，雖然基本能剔除大部分的建筑物等地物點，但是對一些地形復雜含有山地的陡峭地區(qū)時，容易將山頂信息剔除。窗口太小時，容易將一些尺寸較大的建筑物當成地面點保留。如圖6所示。

圖6 建筑物和山地的窗口大小問題

除了窗口大小問題，還有高差閾值無法自動適應地形變化的原因，導致部分非地面點被識別為地面點而產(chǎn)生誤差。比如，有坡度的區(qū)域由于高差閾值相應會增大，會導致斜坡上的房屋植被等非地面點被歸類為地面點。

布料模擬濾波算法對于平坦地形地區(qū)的濾波效果十分良好，但是對于地形較為復雜的區(qū)域(如陡坡、山地等)或者低矮建筑物較多區(qū)域的效果不太好[20]。要么因為平地地區(qū)硬度參數(shù)設置過小或者迭代次數(shù)過多導致把一些低矮建筑物被算成了地面點，要么在陡坡山地地區(qū)因為硬度參數(shù)設置過大或者迭代次數(shù)過少導致山頂或者其他一些地形有變化區(qū)域的地面點被剔除，如圖7所示。

圖7 CSF算法易出現(xiàn)的2類錯誤

2.3.2 高程歸一化

考慮到2種傳統(tǒng)算法各自的缺陷，而布料模擬濾波算法最適應的地形為平坦地形，這種地形下能獲得最好的濾波效果，能夠保留大部分地面點的同時也能有效地過濾非地面點云。因此如何消除高程起伏帶來的影響是關鍵，對點云數(shù)據(jù)進行一個歸一化的操作，使得點云的高程值得到歸一化，來消除崎嶇地形帶來的濾波精度影響。具體方法是先利用原始的點云數(shù)據(jù)粗分類的地面點生成DEM，利用DEM數(shù)據(jù)和原始點云數(shù)據(jù)進行歸一化操作。

考慮到漸進形態(tài)學濾波在對地面點的正確分類方面的有效性，分離出的地面點點云大致反映了目標區(qū)域的地形變化特征，能夠保留地形細節(jié)，因此采用漸進式數(shù)學形態(tài)學濾波所分離出來的地面點構建不規(guī)則三角網(wǎng)，選取點云格網(wǎng)中的最低點為種子點，和最鄰近的2個點構建初始三角形，隨后，各頂點分別尋找最鄰近的一個點組成不規(guī)則三角網(wǎng)(TIN)表面，在三角網(wǎng)中提取柵格格網(wǎng)的值，在三角網(wǎng)中，格網(wǎng)落到三角形頂點時，直接獲得三角形頂點的高度值，如果不在頂點，則該點的高度值通過線性插值方法獲得，在三角形邊的2個頂點上的采用2個頂點的高程內插，在三角形中的則采用各頂點高程內插[21]。

生成的DEM數(shù)據(jù)為TIF格式，將其和原始的點云數(shù)據(jù)結合進行高程歸一化，通過每個點云所處的格網(wǎng)位置找到對應的DEM格網(wǎng)處，點云數(shù)據(jù)的原始高程值和粗DEM對應格網(wǎng)處的高程值的差值就是歸一化的高程值。歸一化之后要保存原始的點云高程值，用于濾波完之后的反歸一化。然后對歸一化之后的點云使用布料模擬濾波，具體的步驟在前面關于布料模擬濾波方法的介紹中已經(jīng)給出，按照原有的步驟進行迭代循環(huán)，分離出地面點和非地面點，濾波結束后賦予原有高程值進行反歸一化，得到結果。技術路線如圖8所示。

圖8 本文技術路線

3 實驗結果及分析

實驗所用的數(shù)據(jù)來自ISPRS官方網(wǎng)站的標準數(shù)據(jù)，選取samp51、samp12、samp52這3塊區(qū)域，均為已經(jīng)經(jīng)過人工編輯分類好的標準數(shù)據(jù)，進行濾波前已經(jīng)進行過粗差剔除，可以直接進行結果分析，一共選擇了3組數(shù)據(jù)。實驗采取C++/PCL+CloudCompare軟件的編譯方式進行，運行環(huán)境為Windows10系統(tǒng)。3個區(qū)域數(shù)據(jù)的基本情況如表1所示。

表1 研究區(qū)域點云基本信息

3．1 實驗結果

samp51數(shù)據(jù)處理過程如圖9所示，其中的內容從左到右、從上到下包括原始的標準點云數(shù)據(jù)、由漸進形態(tài)學濾波結果生成的DEM柵格數(shù)據(jù)、高程歸一化處理后的點云以及最終的濾波結果。而本文方法相比漸進形態(tài)學和普通布料模擬濾波結果的比對如圖10所示，可以看出，本文方法結果使得漸進形態(tài)學濾波中因高差閾值的不確定性或者窗口大小問題而未被剔除的部分非地面點被剔除，同時因普通布料模擬濾波而被錯誤剔除的一些有地形起伏變化的斜坡或者陡峭山地地面點云也得到保留，同理其他2個樣本samp21和samp52的處理結果也類似，samp12結果和比對分別如圖11、圖12所示，samp52結果和比對分別如圖13、圖14所示。

圖9 samp51的濾波處理過程

圖10 samp51和本文算法的比對結果(左形態(tài)學，右布料模擬濾波)

圖11 samp12的濾波處理過程

圖12 samp12和本文算法的比對結果(左布料，右形態(tài)學)

圖13 samp52的濾波處理過程

圖14 samp52和本文算法的比對結果(左布料，右形態(tài)學)

3．2 定量分析

為了驗證算法的精度和可靠性，需要對最終的濾波結果進行一個定量分析。在定量分析中，把誤差分為2類：Ⅰ類誤差和Ⅱ類誤差。Ⅰ類誤差定義是地面點識別為非地面點的誤差，也稱作為拒真誤差，Ⅱ類誤差是非地面點識別為地面點的誤差，也稱作為拒偽誤差。如表2所示，e和f分別表示標準數(shù)據(jù)中的地面點和非地面點的個數(shù)；n為數(shù)據(jù)點的總個數(shù)，而3種誤差的計算公式如下[22]：

(4)

(5)

(6)

表2 濾波誤差定義

將濾波結果和標準數(shù)據(jù)進行比對，本文實驗所選取的參數(shù)在同一樣本下都一致，根據(jù)點云密度范圍確定格網(wǎng)大小，漸進形態(tài)學濾波的參數(shù)如下：坡度值根據(jù)地形變化選取，最小距離閾值為0.5 m，最大距離閾值為3 m，最大窗口值為10 m。布料模擬濾波的參數(shù)如下：距離閾值為0.4 m～0.5 m，迭代次數(shù)為500次，DEM格網(wǎng)大小為2 m或1.5 m。最后得到的結果經(jīng)過計算后各樣本各方法的結果如表3所示。

表3 各區(qū)域濾波實驗誤差統(tǒng)計結果

由表3結果可知，本文方法在進行高程歸一化處理之后的點云再進行布料模擬濾波后結果均比傳統(tǒng)的漸進形態(tài)學濾波和布料模擬濾波要好，相比漸進形態(tài)學濾波，在保證了Ⅰ類誤差穩(wěn)定的同時，降低了Ⅱ類誤差，相比布料模擬濾波，在保證了Ⅱ類誤差穩(wěn)定的同時，降低了Ⅰ類誤差。同時，總誤差也均有降低。

將本文方法和近2年的其他相關改進算法作對比，只選擇其他文中有關本文的3個研究區(qū)域的結果。選取了2020年周欽坤等[23]、2018年陳斐然等[24]改進的數(shù)學形態(tài)學濾波算法、2019年王逸超等[25]一種改進的多級曲面濾波，其結果如表4所示。

表4 濾波結果誤差比較

結果表明，雖然本文算法個別的I類誤差和II類誤差可能相較其他方法更高，但總體而言，總誤差都比其他幾類方法低。

4 結束語

本文對原始點云數(shù)據(jù)進行了一個高程歸一化的操作，解決了傳統(tǒng)的布料模擬濾波和漸進形態(tài)學各自存在的一些問題，同時對濾波結果進行了對比驗證，結果及分析表明：

1)對高程歸一化之后的點云數(shù)據(jù)再進行布料模擬濾波，在地形復雜崎嶇的地方有非常不錯的效果，優(yōu)于傳統(tǒng)的漸進形態(tài)學濾波和布料模擬濾波，也優(yōu)于部分其他的改進算法。

(2)本文方法依賴于生成的DEM精度，因此在處理之前進行點云粗差剔除是必不可少的步驟，本文采用了漸進形態(tài)學濾波作為獲取DEM的前處理過程，也可以嘗試其他算法比如漸進三角網(wǎng)加密。

3)本文需要人工設置DEM格網(wǎng)大小，步驟之間的自動化方面還有待加強。