植被茂密地區(qū)點(diǎn)云雙重濾波方法研究

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

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

引 言

機(jī)載激光雷達(dá)(light detection and ranging,LiDAR)系統(tǒng)是一種集激光測距技術(shù)、全球?qū)Ш叫l(wèi)星定位系統(tǒng)動(dòng)態(tài)定位技術(shù)、慣性導(dǎo)航系統(tǒng)姿態(tài)測量技術(shù)和計(jì)算機(jī)技術(shù)4種技術(shù)于一體的空間測量系統(tǒng)[1]，可獲得高分辨率高精度的空間信息。目前，機(jī)載LiDAR技術(shù)已廣泛應(yīng)用于數(shù)字地面模型獲取、道路提取、電力線提取、林業(yè)調(diào)查、3維城市模型建立、點(diǎn)云分類等地球空間信息學(xué)科的眾多領(lǐng)域[2]。

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

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

1 點(diǎn)云雙重濾波

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

Fig.1 Flow chart of filtering algorithm

1.1 點(diǎn)云偏度平衡-Otsu聯(lián)合粗濾波

1.1.1 基于偏度平衡的單次回波粗濾波 由于原始LiDAR點(diǎn)云中存在噪聲點(diǎn)，首先利用TERRASOLID軟件通過目視檢查剔除部分相對明顯的高位和低位噪聲點(diǎn)，然后結(jié)合強(qiáng)度信息進(jìn)行頻數(shù)統(tǒng)計(jì)分析，剔除強(qiáng)度異常值。剔除噪聲點(diǎn)后，對單次回波的實(shí)驗(yàn)數(shù)據(jù)結(jié)合強(qiáng)度信息利用偏度平衡理論自動(dòng)計(jì)算強(qiáng)度閾值。

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

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

式中，Sk為偏度，μ3為3階中心矩，σ為標(biāo)準(zhǔn)差;(4)若Sk<0，則激光強(qiáng)度閾值I自動(dòng)加1，即I=I+1；若Sk>0，輸出結(jié)果Ith，即點(diǎn)云數(shù)據(jù)的最佳強(qiáng)度閾值;(5)將所有反射強(qiáng)度值大于最佳強(qiáng)度閾值Ith的點(diǎn)云數(shù)據(jù)保留,并參與后續(xù)精濾波。

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

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

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

(2)

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

(3)

將兩式聯(lián)立可得：

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

(4)

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

1.2 點(diǎn)云精濾波

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

具體步驟如下：

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

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

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

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

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

式中，?」為向下取整；xj,k為第j塊第k個(gè)點(diǎn)的x坐標(biāo)；xj,min為第j塊的坐標(biāo)x的最小值；yj,k為第j塊第k個(gè)點(diǎn)的y坐標(biāo)；yj,min為第j塊的坐標(biāo)y的最小值。

(5)遍歷格網(wǎng)。通過索引對每個(gè)格網(wǎng)中點(diǎn)云數(shù)據(jù)的高程值z從小到大進(jìn)行排序，則第1個(gè)數(shù)據(jù)值即為格網(wǎng)中的最小高程值對應(yīng)激光點(diǎn)，所有格網(wǎng)獲取的第一個(gè)數(shù)據(jù)構(gòu)成地面種子點(diǎn)集。

(6)將最后獲取的地面種子點(diǎn)構(gòu)建初始TIN,然后不斷迭代加密實(shí)現(xiàn)點(diǎn)云數(shù)據(jù)精濾波，得到地面點(diǎn)集。

2 實(shí)驗(yàn)數(shù)據(jù)與結(jié)果分析

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

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

Fig.2 Three experimental data sets

2.2 點(diǎn)云粗濾波結(jié)果

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

Fig.3 Schematic diagram of single and first echo intensity

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

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

Fig.4 Experimental data single echo result graph

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

Fig.5 Figure of final echo results of experimental data

基于偏度平衡的單次回波粗濾波和基于Otsu法的多次回波粗濾波分別完成后，得到對應(yīng)單次回波和末次回波的粗濾波結(jié)果并進(jìn)行融合，實(shí)現(xiàn)了偏度平衡-Otsu聯(lián)合粗濾波，如圖6所示。從圖中能明顯看到道路和建筑物被剔除，地形結(jié)構(gòu)得到了較好的保留，為后續(xù)地面種子點(diǎn)的選取提高了可靠性。

Fig.6 Figure of coarse filtering result of point cloud

2.3 點(diǎn)云精濾波結(jié)果

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

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

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

Table 1 Precision evaluation

3類誤差的計(jì)算公式如下：

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

Table 2 Error evaluation table of experimental results

3 結(jié) 論

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