植被茂密地區點云雙重濾波方法研究

王云云，唐菲菲*，王章朋，肖 敏，唐天俊，王銅川

(1.重慶交通大學 土木工程學院，重慶 400074；2.上海寶冶集團有限公司，上海 201908；3.杭州錢塘測繪有限公司，杭州 310000)

引 言

機載激光雷達(light detection and ranging,LiDAR)系統是一種集激光測距技術、全球導航衛星定位系統動態定位技術、慣性導航系統姿態測量技術和計算機技術4種技術于一體的空間測量系統[1]，可獲得高分辨率高精度的空間信息。目前，機載LiDAR技術已廣泛應用于數字地面模型獲取、道路提取、電力線提取、林業調查、3維城市模型建立、點云分類等地球空間信息學科的眾多領域[2]。

盡管機載LiDAR系統硬件已發展較為成熟，然而，機載LiDAR點云數據的濾波算法仍相對比較落后，許多學者對點云濾波算法進行了研究，并提出了多種濾波算法，大概可概括為兩種：第1種是利用點云數據的幾何特性(3維坐標)進行的濾波；第2種是利用點云數據的回波特性(回波強度、多回波信息)進行的濾波[3-6]。經典的濾波算法都是基于幾何特性提出的，如：基于坡度的方法、基于數學形態學的方法、基于曲面擬合的方法以及基于不規則三角網(triangulated irregular network，TIN)算法。近年來，硬件設備快速發展，點云的回波特性得到了學者們的重視。LIU等人利用不同介質的反射率，提出了一種激光強度信息與點云高程信息結合對LiDAR數據進行分類的算法[7]。ZHANG等人提出利用首末次回波的高程差實現了機載LiDAR數據的分類[8]。TANG等人基于點云數據回波信息，將點云以不同分辨率的體素為單位實現點云濾波，在森林地區達到了較好的濾波效果[9-10]。LIN等人利用激光脈沖的多回波特性，預先剔除一部分非地面點的激光腳點，增強了濾波的效果[11]。SUN等人提出了一種基于多回波及Fisher判別的濾波算法，驗證了該方法在陡坡點云濾波中的適應性與有效性[12]。LI等人提出一種基于回波強度的k均值聚類算法，該方法可提高研究區點云濾波精度并減少地面腳點的數據量[13]。

在前人的基礎上，作者提出一種結合幾何特性與回波特性實現點云數據濾波的方法。針對回波特性，知道多次回波會產生不同反射階段，導致激光的能量均有所減弱，就會使單次回波和末次回波之間強度存在差異；且植被密集區域通過多回波分離，發現單次回波的地物類型多樣(道路、建筑、植被)，而末次回波的地物較單一(植被)。因此，本文中對單次回波和末次回波點云數據分別結合不同特征利用不同方法實現自適應粗濾波。另外，不規則三角網漸進加密濾波算法對各種不同地形能取得較好的濾波效果，但其種子點的選取至關重要，所以本文中引入分塊思想，利用局域點云數據密度來選取種子點，使種子點選取更加準確，進而提高不規則三角網漸進加密精濾波算法的結果精度。

1 點云雙重濾波

利用多回波信息進行濾波之前，需根據原始點云的回波次數和第幾次回波數據序列分離出各回波信息，即回波分離。通過回波分離結果分析得出：單次回波包含地面點、大部分植被冠層點以及建筑物等多種地物的回波信號；首次回波則大概率都是植被冠層、樹干等非地面點；末次回波為地表點和植被點的集合；而中間次回波數據量極其少，因此不做分析處理。基于此，本文中只選取單次回波和末次回波的激光腳點參與濾波算法，且針對單次回波和末次回波所包含地物類別的差異性，利用不同的方法分別對單次回波和末次回波進行粗濾波，提出了一種基于偏度平衡-最大類間方差的聯合粗濾波方法；然后，將粗濾波后的單次回波和末次回波數據疊加，參與后續不規則三角網漸進加密精濾波處理，實現點云自適應雙重濾波方法。操作流程如圖1所示。

Fig.1 Flow chart of filtering algorithm

1.1 點云偏度平衡-Otsu聯合粗濾波

1.1.1 基于偏度平衡的單次回波粗濾波 由于原始LiDAR點云中存在噪聲點，首先利用TERRASOLID軟件通過目視檢查剔除部分相對明顯的高位和低位噪聲點，然后結合強度信息進行頻數統計分析，剔除強度異常值。剔除噪聲點后，對單次回波的實驗數據結合強度信息利用偏度平衡理論自動計算強度閾值。

2006年,BARTELS等人最先提出了一種基于偏度平衡的濾波算法[14]，通過引入偏度統計量與峰度統計量的平衡思想，該方法憑借其非監督、無閾值的優點被廣泛使用。根據統計學的中心極限定理，BARTELS提出了兩個假設：(1)自然狀態下，無序離散分布的點云數據中的地面點數據的概率密度函數總是呈正態分布;(2)由于植被、建筑物等非地面點云數據的存在，影響了原始點云數據中的地面點數據的正態分布，導致LiDAR點云數據的概率密度函數從原本近似正態分布呈現出偏態分布[15]。因此，將原始LiDAR點云數據的概率密度函數呈現的偏態分布曲線校正為近似正態分布的過程，就是把點云數據中的非地面點濾除來實現點云數據濾波的過程。本文中基于偏度平衡理論剔除部分非地面點，實現單次回波點云數據的粗濾波。

基于偏度平衡的點云粗濾波方法的具體實現步驟如下：(1)將去噪后的單次回波的點云數據(包括3維坐標xi,yi,zi和強度信息Ii(i=1,…,N)，N為點云總數)單獨標記出來;(2)對點云強度信息進行統計，獲取其最大、最小強度值Imax,Imin，并將最小強度值Imin設為初始強度閾值I;(3)計算單次回波的點云強度偏度值Sk，然后判斷Sk值的正負。計算公式如下：

式中，Sk為偏度，μ3為3階中心矩，σ為標準差;(4)若Sk<0，則激光強度閾值I自動加1，即I=I+1；若Sk>0，輸出結果Ith，即點云數據的最佳強度閾值;(5)將所有反射強度值大于最佳強度閾值Ith的點云數據保留,并參與后續精濾波。

1.1.2 基于最大類間方差的多次回波粗濾波 最大類間方差法是由日本學者OTSU于1979年提出的[16]，是一種圖像分割中確定閾值的最佳算法,簡稱Otsu法。本文中利用該方法獲取最佳閾值的優勢，計算首末次回波的高差閾值，進而實現末次回波點云的粗濾波。

記T為高差閾值，大于高差閾值的點所占比例為ω0,平均高差為μ0；小于高差閾值的點比例為ω1，平均高差為μ1，首末次回波高差的總平均高差為μ，計算類間方差δ，則有[17]：

μ=ω0×μ0+ω1×μ1

(2)

δ=ω0×(μ0-μ)2+ω1×(μ0-μ)2

(3)

將兩式聯立可得：

δ=ω0×ω1×(μ0-μ1)2

(4)

當方差最大時，可以認為獲取的高差為最佳閾值，所以利用最大類間方差法可自適應計算閾值的特性，實現首末次高差的自適應閾值。其具體的算法步驟如下：(1)首先將同一束激光的首次回波和末次回波一一對應，進行標記;(2)計算同束激光對應首末次回波的高程之差，統計其最大、最小值分別為Hmax和Hmin;(3)假設首末次高差進行分類的最佳閾值記為T，則有Hmin

1.2 點云精濾波

經典不規則三角網漸進加密濾波算法是AXELSSON于2000年提出的[18]，該方法對不同場景的適應能力較強、濾波質量較好[19]。通過選取初始地面種子點構建初始不規則三角網，由于初始三角網的精度會影響后續加密過程的判斷，則地面種子點的選取一定程度上會影響濾波的效果[20]。在植被茂密地區，激光點云數據分布很不均勻，因此，本文中引入分塊點云密度，基于局部點云密度設置自適應格網大小，從而獲取可靠地面種子點。

具體步驟如下：

(1)統計粗濾波后獲取的點云數量N和點云數據所占面積S，計算參與精濾波的激光點云分布的整體平均密度ρ=N/S。

(2)點云數據按最大建筑物邊長L進行分塊，記每塊面積為Sj(j=1,2,3,…)，計算分塊后激光點云分布的局部平均密度ρj=Nj/Sj(j=1,2,3,…),其中,Nj為每塊實驗區的點云數量。

(3)點云分布的局部平均密度與整體平均密度ρ進行比較，若ρj>ρ,則計算規則格網大小為gj=int[M÷ρ]-B(j=1,2,3,…)；若ρj<ρ,則計算規則格網大小為gj=[M÷ρ]+B,其中,M表示每個格網中必須包含的最少點數量，B為常數。

(4)點云數據格網分配并建立索引，根據計算得到的虛擬格網大小gj(j=1,2,…)、分塊邊長L以及局部點云最大高差值Hj，可以將局部區域分別確定一個mj×nj×1(j=1,2,…)的虛擬格網，其中,mj和nj的計算方式為：

每塊虛擬格網的編號(Rj,Cj)的計算公式為：

式中，?」為向下取整；xj,k為第j塊第k個點的x坐標；xj,min為第j塊的坐標x的最小值；yj,k為第j塊第k個點的y坐標；yj,min為第j塊的坐標y的最小值。

(5)遍歷格網。通過索引對每個格網中點云數據的高程值z從小到大進行排序，則第1個數據值即為格網中的最小高程值對應激光點，所有格網獲取的第一個數據構成地面種子點集。

(6)將最后獲取的地面種子點構建初始TIN,然后不斷迭代加密實現點云數據精濾波，得到地面點集。

2 實驗數據與結果分析

2.1 實驗數據

本文中選取了3組具有不同地形特征的數據集，如圖2所示。圖2a中的數據集1為平坦區域，包含建筑物、道路、少許植被等地物；圖2b中的數據集2為較平坦區域，但包含大量植被；圖2c中的數據集3為地形起伏區域，包含植被、少許建筑物、道路等地物。3組實驗數據集的單次回波和末次回波點云總數分別為：3882510，3716640，4570200。其中對應單次回波激光腳點數分別為：3721050，3327300，3950280，則末次回波激光腳點數分別為：161460，389340，619920，且每條數據均記錄了3維位置坐標x,y,z,回波號和回波強度信息。

Fig.2 Three experimental data sets

2.2 點云粗濾波結果

2.2.1 基于偏度平衡的單次回波濾波結果 實驗區的單次回波點云數據中主要包含建筑物、道路以及大量植被點等地物，由于不同介質有不同反射率，如：混凝土、瀝青等物質的反射率在20%左右；植被不同種類各不相同，大概在30%～60%范圍；而黏土的反射率則接近于75%，則部分植被與土壤的反射率較為接近[7]，則依據不同地物的反射率存在差異的特性，可利用強度信息濾除一大部分地物點。而在植被覆蓋地區，多回波中首次回波點云數據一般為樹冠、樹干等，且考慮到首次回波沒有激光能量消耗，將其作為單次回波的一個比較對象。圖3為單次回波和首次回波點云強度示意圖。其中條形圖表示單次回波，折線表示首次回波，橫坐標軸為無量綱量強度值，縱坐標為對應強度值點云個數。

Fig.3 Schematic diagram of single and first echo intensity

從單次回波強度分布可以發現，其強度在0～20范圍的數據量極少，而255數據量極大，造成強度變化幅度過大，影響整體回波強度分布規律，則將單次回波強度值為0～20和255的數據作為噪聲點或異常點剔除。另外，圖中單次回波有兩個峰值，與首次回波對應分析可得出，峰值在200左右大概率為地面點，峰值在40左右大概率為地物點，而首次回波在100左右出現數據量急速減少趨勢，可認為強度值在100之前幾乎均為非地面點，對單次回波強度值為100左右的點云進行可視化分析，通過強度數據和可視化判別的綜合分析，最后選取強度值在80～254范圍的單次回波點云數據參與粗濾波。

通過MATLAB平臺編寫該實驗的程序，實現基于偏度平衡的粗濾波方法。將實驗區的單次回波數據序列包含3維坐標x,y,z和強度值(80～254)加載到該濾波方法程序中，作為輸入數據，最后輸出相應的強度閾值分別為160,145,130，結果如圖4所示。從圖中可以看出,利用強度信息能較好地濾除建筑物、道路及大量植被點等地物，且較好保留了地面形態細節，增加了后續種子點選取的準確性。

Fig.4 Experimental data single echo result graph

2.2.2 基于最大類間方差的多次回波濾波結果 實驗數據的首末次回波主要分布在植被覆蓋區域，計算同一激光束對應的首末次回波對應高差值，采用最大類間方差法對實驗數據的首末次回波的高差通過MATLAB編程實現閾值自適應，得到高差閾值分別為6.5m,7.0m,6.0m。將實驗數據中高差大于閾值所對應的末次回波激光腳點提取出來，參與后續濾波處理，其結果如圖5所示。

Fig.5 Figure of final echo results of experimental data

基于偏度平衡的單次回波粗濾波和基于Otsu法的多次回波粗濾波分別完成后，得到對應單次回波和末次回波的粗濾波結果并進行融合，實現了偏度平衡-Otsu聯合粗濾波，如圖6所示。從圖中能明顯看到道路和建筑物被剔除，地形結構得到了較好的保留，為后續地面種子點的選取提高了可靠性。

Fig.6 Figure of coarse filtering result of point cloud

2.3 點云精濾波結果

將粗濾波后的點云數據參與改進的不規則三角網漸進加密精濾波，獲取地面點集。為獲取較高可靠性的地面種子點，首先將實驗區域按照20m×20m進行分塊劃分，然后分別計算分塊后的局部點云密度，將局部點云密度ρi(i=1,2,3,…)與總平均密度ρ進行比較來確定格網大小，其中每個格網中必須包含的最少點數量M和常數B值的設置與平均密度相關，分別分別設置為150,1;100,1;100,2。最后將每個格網中數據按z值進行排序，選取最低點作為地面種子點，進行不規則三角網漸進加密濾波算法實現。圖7為3組數據精濾波得到的地面點俯視圖。從圖中可以看出,該方法較好保留了地形結構，但在植被茂密區域，獲取地面點較小甚至會出現一些空洞，這是以后需要進行改進的一個方向。

Fig.7 Figure of experimental results of point cloud precision filtering

對濾波方法進行結果評價時，僅通過定性分析具有一定片面性，可結合定量分析綜合評定濾波結果。3類誤差是評定點云濾波效果的一個重要指標，將誤差分為第Ⅰ類誤差、第Ⅱ類誤差和總誤差，其中Ⅰ類誤差指地面點誤分為非地面點，Ⅱ類誤差指非地面點誤分為地面點，總誤差為兩類誤差的綜合。表1為精度評定表。表中,e表示參考數據的地面點，即a+b；f表示參考數據的非地面點，即c+d；p表示本文中濾波結果的地面點，即a+c；q表示本文中濾波結果的非地面點，即b+d；h為總數據[21]；且a和d分別為正確濾波的地面點和非地面點，b和c則為錯誤分類的地面點和非地面點。

Table 1 Precision evaluation

3類誤差的計算公式如下：

利用該方法分別對3組實驗數據經過精度進行評定，計算的3類誤差分別如表2所示。從表中可以看出，該濾波方法得到的濾波結果綜合效果較好，數據集1的地形相對平坦且植被覆蓋較低，其3類誤差都相對較小；數據集2則植被覆蓋率極高，相比其它兩組數據集的3類誤差相對較大，但第Ⅱ類誤差即誤分為地面點的誤差相對來說也是較低的，說明獲取的地面點準確度相對比較高；數據集3的地形有一定坡度且植被覆蓋率較高，其結果精度相對較高，也較好保留了地形細節。因此，本文中提出的雙重濾波方法不僅在平坦且植被覆蓋較少的地區取得較好結果，而且在地形起伏且茂密植被覆蓋區域仍能較好保留地形細節，但植被覆蓋率過高時，可考慮與其它多源數據進行融合分析。

Table 2 Error evaluation table of experimental results

3 結 論

本文中的方法對單次回波和末次回波的激光腳點進行粗濾波，實現閾值自動化，且降低了后續種子點選取無法自動剔除建筑物和植被等非地面點的可能性，大大提高了精濾波運算效率和濾波精度，但該方法只適用于植被高度分布相對均勻的區域，高差閾值的確定方法需要進一步研究。另外，針對植被茂密地區濾波后出現的稀疏點云或者空洞區域如何處理仍需深入研究。