程鐵洪廖永福胡冬清

（1.江西省電力設(shè)計(jì)院 江西南昌 330096；2.江西省交通設(shè)計(jì)研究院 江西南昌 330209）

1 概述

隨著機(jī)載激光雷達(dá) （Light Detection and Ranging，LiDAR）技術(shù)的發(fā)展，它在各行各業(yè)得到了廣泛的應(yīng)用，尤其是機(jī)載激光雷達(dá)能快速準(zhǔn)確的獲取地面物體的三維坐標(biāo)，使得生成大規(guī)模、高精度的DEM已經(jīng)越來越方便。但原始點(diǎn)云中不僅含有地面點(diǎn)，還有較多的非地面點(diǎn)。因此，通過機(jī)載激光雷達(dá)技術(shù)獲取地面點(diǎn)的三維坐標(biāo)，進(jìn)而生成DEM，首先要解決的就是濾波問題，LiDAR點(diǎn)云數(shù)據(jù)濾波就是從原始的LiDAR點(diǎn)云中濾掉非地面點(diǎn)，保留有效的地形信息。目前主要的濾波算法包括Lindenberger［1］和 Kilian［2］等人提出的數(shù)學(xué)形態(tài)學(xué)濾波算法、Axelsson［3］提出逐漸加密三角網(wǎng)(TIN)濾波算法以及Kraus和 Pfeifer［4］提出的線性預(yù)測(cè)濾波算法等。

大部分的機(jī)載LiDAR點(diǎn)云濾波算法的基本原理是基于鄰近激光腳點(diǎn)高程突變的原理，基于該原理有以下幾個(gè)基本假設(shè)：（1）一定范圍內(nèi)的局部區(qū)域高程最低點(diǎn)為地面點(diǎn)；（2）掃描的區(qū)域相對(duì)比較平緩，沒有劇烈的地形起伏［5-6］；（3）地形曲面由一片片連續(xù)光滑的空間曲面所構(gòu)成。各種濾波算法都有一定的局限性，而線性預(yù)測(cè)濾波算法對(duì)不同地形的適應(yīng)性較好，處理結(jié)果精度較高，因此被廣泛的使用。

傳統(tǒng)的線性預(yù)測(cè)濾波算法仍存在以下不足：（1）噪聲點(diǎn)的影響，噪聲點(diǎn)對(duì)后期地形趨勢(shì)面的擬合效果和濾波參數(shù)設(shè)置都有較大影響，原算法未考慮噪聲點(diǎn)的影響；（2）數(shù)據(jù)分塊大小固定，原算法未考慮不同地形條件應(yīng)合理分塊的情況；（3）濾波循環(huán)終止條件不夠明確；（4）濾波參數(shù)計(jì)算復(fù)雜，當(dāng)濾波效果不佳時(shí)經(jīng)常需要對(duì)參數(shù)進(jìn)行反復(fù)的調(diào)試，人工干預(yù)性較大；（5）在森林茂密地區(qū)，當(dāng)?shù)孛纥c(diǎn)較少時(shí)，地形趨勢(shì)面擬合效果不好，很難獲得林區(qū)真實(shí)的地形情況。

針對(duì)上述不足，對(duì)傳統(tǒng)線性預(yù)測(cè)在以下幾個(gè)方面進(jìn)行了改進(jìn)：（1）噪聲點(diǎn)剔除，通過分塊內(nèi)高程分布直方圖對(duì)點(diǎn)云進(jìn)行了有效的去噪；（2）數(shù)據(jù)格網(wǎng)化分塊，根據(jù)地形條件和點(diǎn)云密度不同對(duì)點(diǎn)云進(jìn)行分塊處理，能有效提高濾波效率；（3）迭代終止條件判斷，通過引入偏度的概念，計(jì)算迭代過程中點(diǎn)云殘差的偏度sk值，可以有效的判定迭代次數(shù)；（4）濾波參數(shù)設(shè)置，通過分析殘差分布直方圖，確定濾波過程中的重要參數(shù)g和w，提高了算法對(duì)地形的適應(yīng)性。

2 算法原理及技術(shù)流程

2.1 傳統(tǒng)的線性預(yù)測(cè)濾波算法

傳統(tǒng)的線性預(yù)測(cè)算法在進(jìn)行濾波時(shí)，首先將原始數(shù)據(jù)劃分為若干小塊，通過對(duì)分塊內(nèi)點(diǎn)云擬合參考表面，這個(gè)參考表面一般為二次曲面［7］，其解析式為：

式中，i為點(diǎn)號(hào)，x，y為平面坐標(biāo)，z為擬合后的點(diǎn)高程，a0至a5為擬合曲面參數(shù)。

然后用原始數(shù)據(jù)各點(diǎn)的高程值減去對(duì)應(yīng)參考表面的擬合值，得到該點(diǎn)的擬合殘差，根據(jù)殘差確定該點(diǎn)在下一次迭代計(jì)算中的權(quán)值，權(quán)函數(shù)關(guān)系式為：

式中，pi表示第i個(gè)點(diǎn)的權(quán)值，vi表示第i個(gè)點(diǎn)的殘差，參數(shù)a，b的值決定了權(quán)函數(shù)的陡峭程度，通常依靠經(jīng)驗(yàn)a，b可分別取1和4，w為某分塊區(qū)域內(nèi)地形的最大起伏。在一次迭代過程中，依據(jù)權(quán)函數(shù)關(guān)系式和殘差對(duì)點(diǎn)云進(jìn)行分類時(shí)，引入了濾波重要參數(shù)偏移量g和w，用于判斷點(diǎn)云的屬性，權(quán)函數(shù)示意圖如下圖1所示：

圖1 權(quán)函數(shù)示意圖

參數(shù)g的計(jì)算傳統(tǒng)線性預(yù)測(cè)算法主要通過殘差分布直方圖計(jì)算，方法主要有：（1）用事先預(yù)期的地面點(diǎn)測(cè)量的精度來計(jì)算，從殘差分布直方圖的原點(diǎn)（殘差為0）往左移，直到計(jì)算出的標(biāo)準(zhǔn)差與預(yù)先預(yù)期的地面點(diǎn)精度σT相同時(shí)，這個(gè)移動(dòng)量就是的g值，如下圖 2（a）部分所示；（2）從殘差累計(jì)曲線圖上計(jì)算出向左移動(dòng)每個(gè)位置的標(biāo)準(zhǔn)差，然后繪制累計(jì)標(biāo)準(zhǔn)差曲線圖，如下圖2（b）部分所示，找出標(biāo)準(zhǔn)差達(dá)到極小值時(shí)，對(duì)應(yīng)的向左移動(dòng)的距離為g的值（3）根據(jù)激光雷達(dá)的穿透率，g的值就是殘差直方圖移動(dòng)到地面點(diǎn)有一半處的值。

圖2 計(jì)算g值所用方法示意圖

2.2 改進(jìn)后的算法流程

改進(jìn)后的線性預(yù)測(cè)濾波算法流程如下圖3所示，主要在噪聲點(diǎn)剔除、數(shù)據(jù)格網(wǎng)化分塊、迭代終止條件判斷以及重要濾波參數(shù)設(shè)置上進(jìn)行了改進(jìn)。

圖3 改進(jìn)后的線性預(yù)測(cè)濾波算法流程圖

2.3 改進(jìn)算法的關(guān)鍵技術(shù)

針對(duì)傳統(tǒng)線性預(yù)測(cè)濾波算法存在的不足，做出以下改進(jìn)：

（1）點(diǎn)云剔噪。根據(jù)噪聲點(diǎn)和非噪聲點(diǎn)在空間分布的特點(diǎn)，本文研究采用基于高程分布特點(diǎn)的點(diǎn)云去噪方法。該算法的基本思想是基于噪聲點(diǎn)相對(duì)于局部鄰近點(diǎn)在高程上存在顯著的突變現(xiàn)象，非噪聲點(diǎn)在高程分布上處于相對(duì)連續(xù)的狀態(tài)，這也是目前大多數(shù)機(jī)載LiDAR點(diǎn)云去噪算法的基本假設(shè)。該算法需要對(duì)數(shù)據(jù)集進(jìn)行分塊，在每個(gè)分塊內(nèi)繪制所有點(diǎn)的高程分布直方圖，如下圖4所示，高程的間隔盡量要小，如1米較為適宜，從圖中可以看出，非噪聲點(diǎn)在高程上呈現(xiàn)連續(xù)分布狀態(tài)，噪聲點(diǎn)則呈現(xiàn)離散分布且與非噪聲點(diǎn)分隔開，通過高程分布直方圖可以將噪聲點(diǎn)和非噪聲點(diǎn)在空間的分布特點(diǎn)很好的展示出來，從而將噪聲點(diǎn)剔除。

圖4 sample11區(qū)某分塊內(nèi)點(diǎn)云高程分布統(tǒng)計(jì)直方圖

（2）點(diǎn)云數(shù)據(jù)格網(wǎng)化。在進(jìn)行數(shù)據(jù)格網(wǎng)化分塊時(shí)，主要解決兩個(gè)核心的問題，一是格網(wǎng)大小，另一個(gè)是建立區(qū)塊索引機(jī)制。濾波時(shí)，一個(gè)格網(wǎng)就是一個(gè)濾波窗口，在考慮格網(wǎng)寬度時(shí)，主要是考慮窗口的大小不能小于窗口內(nèi)最大地物面積，若濾波窗口小于最大地物區(qū)域，該窗口的濾波結(jié)果就會(huì)出現(xiàn)系統(tǒng)性誤差，在城區(qū)地勢(shì)較為平坦的地形，格網(wǎng)的寬度一般可設(shè)置為 40×40m［7］，在大型建筑物密集的區(qū)塊，設(shè)置為40～100m不等。依據(jù)實(shí)際情況做出調(diào)整，在地形陡峭的山區(qū)，格網(wǎng)的寬度應(yīng)適當(dāng)減小。建立區(qū)塊索引機(jī)制是對(duì)分塊后點(diǎn)云進(jìn)行濾波時(shí)的快速定位，并且將所在該區(qū)塊的點(diǎn)云集中處理，通常根據(jù)分塊后各點(diǎn)所在的行列號(hào)來進(jìn)行索引。

具體來說，在數(shù)據(jù)格網(wǎng)化分塊時(shí)，具體流程如下：

（1）讀取數(shù)據(jù)預(yù)處理后的點(diǎn)云數(shù)據(jù)；

（2）索引整個(gè)區(qū)域內(nèi)平面坐標(biāo)X和Y的最大最小值 Xmax、Xmin 以及 Ymax、Ymin；

（4）再根據(jù)各點(diǎn)的平面坐標(biāo)將點(diǎn)云分到各個(gè)格網(wǎng)中并獲得它們所在的行列號(hào)，分塊的同時(shí)，計(jì)算落入各個(gè)格網(wǎng)的點(diǎn)云數(shù)。

格網(wǎng)索引信息獲取由公式（3）計(jì)算得到，具體如下：

公式（3）中，row和column分別表示整個(gè)區(qū)域格網(wǎng)的行列數(shù)，c表示該區(qū)域內(nèi)格網(wǎng)的寬度，Xmax、Xmin以及Ymax、Ymin分別表示整個(gè)區(qū)域內(nèi)平面坐標(biāo)的最大最小值，floor表示對(duì)計(jì)算結(jié)果取整數(shù)，xp、yp分別表示格網(wǎng)內(nèi)任意點(diǎn)的平面坐標(biāo)，i、j則分別表示該任意點(diǎn)所在格網(wǎng)的行列號(hào)。點(diǎn)云數(shù)據(jù)格網(wǎng)化示意圖如下圖5所示：

圖5 點(diǎn)云數(shù)據(jù)格網(wǎng)化示意圖

（3）迭代終止條件判斷。引入偏度的概念用于判斷迭代終止條件時(shí)，由數(shù)理統(tǒng)計(jì)中的中心極限定理可知，只含地面點(diǎn)的殘差分布的概率密度函數(shù)會(huì)服從正態(tài)分布［8］，點(diǎn)云中非地面點(diǎn)的存在干擾了原有的正態(tài)分布，使得點(diǎn)云的整體分布呈現(xiàn)某種偏態(tài)分布。因此，濾波的過程就是要剔除非地面點(diǎn)，只保留地面點(diǎn)，也就是使受到非地面點(diǎn)干擾的正偏態(tài)分布重新變?yōu)檎龖B(tài)分布的過程，濾波直到sk=0或接近0即可結(jié)束。在實(shí)際的判斷過程中發(fā)現(xiàn)，由于非地面點(diǎn)無法完全剔除，濾波后殘差的概論密度函數(shù)并不會(huì)完全達(dá)到正態(tài)分布，也就是說偏度值很難達(dá)到0，如下圖6所示，表示某平坦地區(qū)濾波結(jié)束后殘差分布直方圖和概率密度函數(shù)曲線，計(jì)算得偏度sk=0.43，概率密度函數(shù)曲線近似的服從正態(tài)分布，此時(shí)，該區(qū)域?yàn)V波已達(dá)到比較好的結(jié)果。因此，在判斷時(shí)，根據(jù)實(shí)際地形考慮濾波難度時(shí)對(duì)迭代終止時(shí)的偏度值設(shè)置了一個(gè)閾值 sk閾值（sk閾值>0），若某次迭代結(jié)束殘差分布的偏度的絕對(duì)值|sk|小于偏度閾值sk閾值，即|sk|

①平坦地區(qū)：0=

②地形起伏大山區(qū)：0=

圖6 濾波后地面點(diǎn)殘差頻數(shù)分布直方圖

改進(jìn)后的迭代終止判斷條件，無需考慮迭代次數(shù)，無需判斷兩次迭代間地形趨勢(shì)面是否變化，避免了較大的人工干預(yù)。

（4）濾波參數(shù)g和w的確定。利用殘差分布直方圖來確定濾波參數(shù)，首先，要分析殘差分布直方圖的特點(diǎn)，本文從不同地形條件出發(fā)，總結(jié)了不同地形條件下殘差分布直方圖的共同特點(diǎn)，如下圖7所示，主要特點(diǎn)表現(xiàn)為以下幾個(gè)方面：

①從整體來看，殘差分布直方圖不服從正態(tài)分布，而是向非地面點(diǎn)（正殘差部分）偏移，整體呈正偏態(tài)分布；

1.2.3 實(shí)驗(yàn)儀器和試劑 美國Waters生產(chǎn)的ACQUITYTM型串聯(lián)質(zhì)譜儀、PE公司生產(chǎn)的非衍生化多種氨基酸、肉堿和琥珀酰丙酮測(cè)定試劑盒(以后簡(jiǎn)稱試劑盒)、艾本森MicBio-IV型孵育震蕩器、Vortex Mixer XW-80A混勻器以及打孔鉗(打孔柱直徑3 mm)。

②地面點(diǎn)（負(fù)殘差部分）的殘差分布比較集中，呈近似正態(tài)分布，非地面點(diǎn)分布較零散，且殘差相差較大；

③在復(fù)雜地形條件下，地面點(diǎn)殘差分布呈多峰現(xiàn)象。

圖7 不同地形條件下的殘差分布直方圖（首次擬合）

采用原算法中提到的第三種方法，即：g的值就是殘差直方圖原點(diǎn)出發(fā)，向左移動(dòng)到地面點(diǎn)的殘差占整個(gè)殘差一半位置時(shí)的值，移動(dòng)的距離即為g的值。在通過殘差分布直方圖分析w的取值時(shí)，首先，w是用來區(qū)分哪些點(diǎn)為非地面點(diǎn)的外閾值，從殘差分布直方圖表現(xiàn)為正態(tài)分布曲線右半部分趨于0或殘差發(fā)生新的突變處的地方，如下圖8所示。

圖8 通過殘差分布直方圖分析g和w取值方法

3 實(shí)驗(yàn)與分析

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

本文實(shí)驗(yàn)所用數(shù)據(jù)均來自國際攝影測(cè)量和遙感學(xué)會(huì)（ISPRS）在線發(fā)布的參考數(shù)據(jù)作為實(shí)驗(yàn)數(shù)據(jù)，該數(shù)據(jù)對(duì)每個(gè)點(diǎn)都做了人工分類，標(biāo)定地面點(diǎn)和非地面點(diǎn)。本文實(shí)驗(yàn)選取了其中4處比較有代表性的地形進(jìn)行實(shí)驗(yàn)，Sample11（陡峭山坡）、Sample12（密集建筑物）、Sample51（茂密的林區(qū)）以及 Sample71（大型構(gòu)筑物）。

3.2 濾波結(jié)果分析

（1）定性分析

下圖9至圖12是4塊實(shí)驗(yàn)區(qū)域?yàn)V波前后點(diǎn)云生成的DEM及濾波后地面點(diǎn)的高程頻數(shù)直方圖，其中（a）表示濾波前點(diǎn)云生成的DEM，（b）表示人工分類地面點(diǎn)生成的參考DEM，（c）表示濾波后點(diǎn)云生成的DEM，（d）表示濾波后地面點(diǎn)生成的高程頻數(shù)直方圖。

圖9 Sample11濾波前后點(diǎn)云生成的DEM及濾波后地面點(diǎn)高程頻數(shù)直方圖（sk=2.78）

圖10 Sample12濾波前后點(diǎn)云生成的DEM及濾波后地面點(diǎn)高程頻數(shù)直方圖（sk=1.34）

圖11 Sample51濾波前后點(diǎn)云生成的DEM及濾波后地面點(diǎn)高程頻數(shù)直方圖（sk=1.25）

圖12 Sample71濾波前后點(diǎn)云生成的DEM及濾波后地面點(diǎn)高程頻數(shù)直方圖（sk=1.92）

通過對(duì)實(shí)驗(yàn)區(qū)定性分析可以看出：在地形較陡峭的山區(qū)，Sample11濾波效果不佳；在地形較平坦的城區(qū)和林區(qū)，Sample12和Sample51濾波精度較高；在有大型構(gòu)筑物的地形條件下，Sample71存在地形過分割現(xiàn)象，但總體保留了地形的趨勢(shì)面，說明改進(jìn)后的算法具有較好的適應(yīng)性，并且從濾波后殘差的偏度值也能反映濾波總結(jié)較好。

（2）定量分析

對(duì)四個(gè)實(shí)驗(yàn)區(qū)域 Sample11、Sample12、Sample51以及Sample71進(jìn)行濾波算法的定量分析。用I類誤差（將地面點(diǎn)漏分為地物點(diǎn)的錯(cuò)誤率），II類誤差（將地物點(diǎn)錯(cuò)分為地面點(diǎn)的錯(cuò)誤率總誤差）以及總誤差（反應(yīng)將點(diǎn)云錯(cuò)分漏分加權(quán)后的總誤差）來評(píng)定濾波的準(zhǔn)確性，算法改進(jìn)前后三類誤差對(duì)比圖如下圖13所示。

通過對(duì)實(shí)驗(yàn)區(qū)定量分析可以看出：改進(jìn)后的算法I類誤差總體減小，尤其是在地勢(shì)坡度大，起伏度大的地方，如Sample11和Sample71區(qū)域，在地勢(shì)較平坦或者起伏不大的區(qū)域，如Sample12和Sample51區(qū)域，I類誤差和原算法相當(dāng)；雖然造成II類誤差有所增大，但I(xiàn)I類誤差的點(diǎn)在可視化條件下通過人工干預(yù)較容易剔除；從總誤差來看，改進(jìn)后的算法總誤差明顯降低，尤其是在坡度和地勢(shì)起伏較大的地方，總誤差的大小反映濾波算法的可行性。因此，從定量分析角度來看，改進(jìn)后的算法在部分區(qū)域優(yōu)于原算法，一定程度上提高了濾波精度，較好的保留了地形的趨勢(shì)面。

圖13 算法改進(jìn)前后三類誤差錯(cuò)誤率比較

4 結(jié)論

通過對(duì)改進(jìn)后的濾波算法進(jìn)行了定性和定量分析，從分析的結(jié)果來看，改進(jìn)后的算法有以下特點(diǎn)和優(yōu)勢(shì)：

（1）在地形平坦或地勢(shì)起伏不大的區(qū)域，改進(jìn)后的算法和原算法濾波效果相當(dāng)，且都達(dá)到較高的精度。

（2）在坡度較大或地勢(shì)起伏劇烈的區(qū)域，原算法的濾波效果顯著降低，改進(jìn)后的濾波算法更加注重對(duì)地面點(diǎn)的保留，雖然II類誤差有所增大，但總誤差卻有所下降，且更好的保留了地形的趨勢(shì)面，使得濾波后生成的DEM更加符合真實(shí)地形，顯示出改進(jìn)后的算法精度更高，增強(qiáng)了算法的適應(yīng)性。

（3）改進(jìn)的算法增加了濾波后點(diǎn)云高程分布圖的偏度計(jì)算來評(píng)估濾波的質(zhì)量，這在一定程度上能輔助對(duì)濾波結(jié)果的評(píng)估，尤其是在原始點(diǎn)云無人工分類的情況下，這種方法切實(shí)可行。