張 靖,張曉君,江萬(wàn)壽,王建超,郭大海

(1.武漢大學(xué)測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,武漢 430079;2.中國(guó)國(guó)土資源航空物探遙感中心,北京 100083)

一種改進(jìn)的線(xiàn)性預(yù)測(cè)濾波算法

張 靖1,張曉君1,江萬(wàn)壽1,王建超2,郭大海2

(1.武漢大學(xué)測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,武漢 430079;2.中國(guó)國(guó)土資源航空物探遙感中心,北京 100083)

對(duì)傳統(tǒng)的線(xiàn)性預(yù)測(cè)算法進(jìn)行改進(jìn),通過(guò)粗差剔除、初始點(diǎn)選取、特殊地形分析及光滑度檢驗(yàn)等步驟提高了濾波精度。利用 ISPRS發(fā)布的 6組標(biāo)準(zhǔn)測(cè)試數(shù)據(jù)與原方法進(jìn)行實(shí)驗(yàn)比較,結(jié)果證明,改進(jìn)后的算法提高了濾波總精度以及地面點(diǎn)分類(lèi)精度。

L iDAR;濾波;線(xiàn)性預(yù)測(cè);粗差;分塊

0 引言

機(jī)載激光雷達(dá) (LightDetection And Ranging,Li-DAR)是一種快速獲取高精度三維信息的新型遙感技術(shù)。LiDAR點(diǎn)云的濾波就是要從原始 LiDAR數(shù)據(jù)中濾掉非地面點(diǎn),保留有效的地形信息。目前國(guó)際上已提出了多種濾波算法[1],如 Linderberger[2]、Killian等[3]和 Zhang等[4]提出的數(shù)學(xué)形態(tài)學(xué)方法,Axelsson[5]提出的漸進(jìn)式 TIN濾波算法以及 Kraus等[6,7]提出的線(xiàn)性預(yù)測(cè)濾波算法等。

現(xiàn)有濾波方法的基本原理都是基于局部地形的連續(xù)性變化這一假設(shè),算法設(shè)計(jì)的關(guān)鍵在于如何快速準(zhǔn)確地檢測(cè)出局部的高程突變,這需要解決兩個(gè)關(guān)鍵問(wèn)題:①如何確定激光點(diǎn)之間的鄰域關(guān)系;②選擇什么標(biāo)準(zhǔn)作為高程變化的判別條件。由于線(xiàn)性預(yù)測(cè)算法不需要確定激光點(diǎn)之間的嚴(yán)格位置關(guān)系,也不需要人為設(shè)定精確閾值,算法的適應(yīng)性好,處理精度較高,因此被廣泛地采用[8-10],并成為一種標(biāo)準(zhǔn)的L iDAR濾波方法。

但是,傳統(tǒng)的線(xiàn)性預(yù)測(cè)算法也存在著一些缺陷:①分塊大小固定,未考慮分塊方式的科學(xué)性、地形形態(tài)的適應(yīng)性及相鄰分塊的光滑程度;②算法易受低點(diǎn)誤差影響。假設(shè)高程較低的點(diǎn)為地面點(diǎn),當(dāng)數(shù)據(jù)中存在低于地面的奇異點(diǎn)時(shí),會(huì)導(dǎo)致計(jì)算結(jié)果出錯(cuò);③算法處理過(guò)程中需進(jìn)行多次迭代。迭代過(guò)程中用整個(gè)數(shù)據(jù)集擬合趨勢(shì)面,計(jì)算效率較低;④在森林和城鎮(zhèn)地區(qū)地面點(diǎn)較少時(shí),很難估計(jì)出準(zhǔn)確的地面[8]。

針對(duì)上述問(wèn)題,本文提出了一種改進(jìn)的線(xiàn)性預(yù)測(cè)濾波方法。具體改進(jìn)如下:①對(duì)數(shù)據(jù)進(jìn)行分塊時(shí),依照點(diǎn)密度的大小設(shè)置不同的分塊大小,并通過(guò)高程信息的豐富程度劃分不同的層次;②在濾波前通過(guò)對(duì)分塊內(nèi)直方圖進(jìn)行分析,去除奇異點(diǎn);③初次擬合時(shí),選取高程分布在 [15%,65%]范圍內(nèi)的點(diǎn)作初始值,再逐漸加點(diǎn),以減少擬合的計(jì)算量;④通過(guò)分析特殊地形的高程直方圖,對(duì)僅含有房屋點(diǎn)和少量地面點(diǎn)的分塊進(jìn)行擬合初值設(shè)定,提高了算法的適用性。實(shí)驗(yàn)證明,改進(jìn)的線(xiàn)性濾波算法與原算法比較,具有更強(qiáng)的適應(yīng)性和更高的濾波精度。

1 算法描述

1.1 傳統(tǒng)線(xiàn)性預(yù)測(cè)算法

線(xiàn)性預(yù)測(cè)濾波是在濾波的同時(shí)內(nèi)插DEM。對(duì)于給定的一系列離散點(diǎn)集,為了更好地進(jìn)行線(xiàn)性估計(jì),首先將原始數(shù)據(jù)劃分為小塊,用一個(gè)移動(dòng)曲面來(lái)逼近局部區(qū)域的趨勢(shì)面,一般設(shè)為二次曲面,其解析式為

式中,i為分塊號(hào);a1i等為第 i塊曲面的參數(shù);X,Y為平面坐標(biāo)。

然后,用原始數(shù)據(jù)中每一點(diǎn)的高程值減去這點(diǎn)趨勢(shì)面的高程值,即為擬合的殘差。利用擬合殘差確定該點(diǎn)在下一次曲面擬合中的權(quán)重,權(quán)函數(shù)定義為

式中,vi為第 i點(diǎn)的擬合殘差;ρi為第 i點(diǎn)的權(quán)值;g為偏移值;w為地形最大起伏;a,b決定了權(quán)函數(shù)的陡峭程度,其計(jì)算方法為

式中,h為權(quán)值函數(shù)半寬值;s為半寬處的斜率;h和 s的含義如圖 1所示。在一般情況下,可以取a=1,b=4。

圖 1 權(quán)函數(shù)示意圖Fig.1 Diagram of the weight function

參數(shù) g可以根據(jù)殘差直方圖統(tǒng)計(jì)得到,主要有3種方法:①用預(yù)期的地面點(diǎn)精度計(jì)算,將殘差直方圖不斷向左移,直至計(jì)算出的地面點(diǎn)精度與預(yù)期精度相同時(shí),向左移動(dòng)的距離即為 g(但當(dāng)存在大的粗差時(shí),該方法就會(huì)失效);②計(jì)算出向左移動(dòng)每一個(gè)位置后點(diǎn)的標(biāo)準(zhǔn)差,找出最小的標(biāo)準(zhǔn)差,對(duì)應(yīng)的向左移動(dòng)的距離為 g(圖 2);③通過(guò)估計(jì)的地面點(diǎn)的比例來(lái)計(jì)算。如果初始估計(jì)地面點(diǎn)比例為 40%,殘差直方圖要移動(dòng)到有 20%地面點(diǎn)的地方[6]。

圖 2 移動(dòng)距離 g的計(jì)算方法示意圖Fig.2 Calculate moving distanceg

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

圖 3為改進(jìn)的線(xiàn)性預(yù)測(cè)濾波的算法流程圖,其中實(shí)線(xiàn)框標(biāo)明的部分表示傳統(tǒng)的線(xiàn)性預(yù)測(cè)濾波流程,虛線(xiàn)框部分為改進(jìn)的部分。

圖 3 改進(jìn)的線(xiàn)性濾波流程Fig.3 Flow chart of progressive linear prediction filtering

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

針對(duì)上述線(xiàn)性預(yù)測(cè)濾波算法的幾個(gè)缺點(diǎn),本文做出如下改進(jìn):

(1)粗差剔除。一般的濾波算法都認(rèn)為在局部范圍內(nèi)最低點(diǎn)為地面點(diǎn),而機(jī)載 L iDAR在飛行過(guò)程中都會(huì)出現(xiàn)粗差點(diǎn) (即極低點(diǎn)或極高點(diǎn)),許多濾波算法都沒(méi)有去除極低點(diǎn)的能力。本文從計(jì)算高程直方圖的角度簡(jiǎn)單去除粗差。無(wú)論是極低點(diǎn)還是極高點(diǎn),都表現(xiàn)為其高程與其他點(diǎn)的高程差異較大,且個(gè)數(shù)較少。在計(jì)算分塊內(nèi)高程直方圖時(shí),若把份數(shù)劃分得較多,粗差點(diǎn)與正常點(diǎn)之間就會(huì)存在大量的空白頻數(shù),因此可通過(guò)去除空白頻數(shù)之前的點(diǎn),達(dá)到去除低點(diǎn)誤差的目的。

(2)合理分塊。線(xiàn)性預(yù)測(cè)濾波是用局部擬合的方法逐漸逼近地面的,因此分塊不可太大;但為了更好地濾除建筑物,又要求分塊不可太小。本文依據(jù)經(jīng)驗(yàn)估算出最大建筑物的尺寸,一般設(shè)置為40m×40m,并根據(jù)點(diǎn)密度的不同自動(dòng)調(diào)整分塊大小。若窗口內(nèi)點(diǎn)云的密度超過(guò)測(cè)區(qū)的平均點(diǎn)云密度,則將窗口分為 4個(gè)小窗口。這樣,在保證分塊大小相似的同時(shí),確保了每個(gè)分塊內(nèi)有近似數(shù)量的激光點(diǎn),擬合時(shí)可以得到更穩(wěn)定的結(jié)果。

(3)初始點(diǎn)選取。Kraus等[7]提出的迭代線(xiàn)性最小二乘內(nèi)插法,其計(jì)算始終是面向整個(gè)數(shù)據(jù)集,運(yùn)算速度較慢。后來(lái),有人提出首先計(jì)算出整個(gè)區(qū)域內(nèi)所有點(diǎn)的高程平均值 Zi,設(shè)一個(gè)水平面 Z=Zi和一個(gè)高程閾值 a,將位于 Z=Zi-a與 Z=Zi+a之間的點(diǎn)作為初始點(diǎn)進(jìn)行擬合,再加入符合條件的點(diǎn),不斷逼近 (其中 a的域值選取多為分塊內(nèi)最大高程值與最小高程值差值的 1/3～1/4)。這種方法減少了初始擬合點(diǎn)數(shù),提高了運(yùn)算速度,但在有些情況下這一范圍的點(diǎn)太少,會(huì)影響平面擬合的精度?；诖?本文利用高程直方圖統(tǒng)計(jì)的方法來(lái)選擇初始點(diǎn)。一般選取高程數(shù)目分布在[15%,65%]之間的點(diǎn)數(shù)進(jìn)行初步擬合,這樣既可以提高運(yùn)算速度,也可以選到足夠數(shù)量的可靠點(diǎn) (如圖 4所示,可以選取高程范圍在[334,356]范圍內(nèi)的點(diǎn)作為初始點(diǎn))。

圖 4 分塊內(nèi)高程直方圖Fig.4 Histogram of elevation

(4)特殊地形分析。線(xiàn)性迭代最小二乘內(nèi)插法是通過(guò)逐步擬合平面來(lái)逼近地面的,但只有分塊內(nèi)地面點(diǎn)相對(duì)較多時(shí),擬合出來(lái)的平面才能代表正確的地面。有時(shí)房屋面積遠(yuǎn)大于分塊面積,由于分塊內(nèi)地面點(diǎn)數(shù)太少,擬合出的最終平面是介于地面點(diǎn)和房屋頂面中間的平面,使得該方法失效。但從高程直方圖的角度來(lái)看,分塊的選擇一般需滿(mǎn)足以下條件:①高程有兩個(gè)峰值;②峰值間高程相差較大;③峰值相鄰的高程信息較為簡(jiǎn)單;④第一個(gè)峰值的點(diǎn)數(shù)較少。因此可以先判斷分塊是否滿(mǎn)足特殊地形條件,若滿(mǎn)足則選取第一個(gè)峰值附近的點(diǎn)為地面初始點(diǎn)。這樣擬合出來(lái)的平面為真正地面,可以得到較為準(zhǔn)確的濾波結(jié)果。

(5)光滑度檢驗(yàn)。最小二乘內(nèi)插法利用局部擬合平面的方法來(lái)逐步濾出地面點(diǎn)。由于地表被劃分為一個(gè)個(gè)的小塊,因此不僅小塊內(nèi)部要滿(mǎn)足光滑、連續(xù)的條件,分塊間也要滿(mǎn)足光滑的條件。本文先對(duì)每一小塊進(jìn)行擬合濾波,得到該塊的地面表達(dá)式后(如圖 4中的陰影塊),檢驗(yàn)它與 4個(gè)鄰域小塊 (圖 4中的 1,2,3,4)的光滑程度,可以按照不同的要求選擇 0階光滑 (即 f(x,y)=g(x,y))或一階光滑 (即f′(x,y)=g′(x,y))等判斷標(biāo)準(zhǔn)。若不滿(mǎn)足光滑條件,可聯(lián)合陰影部分與該塊共同擬合平面,以提高濾波準(zhǔn)確度。

圖 5 相鄰分塊示意圖Fig.5 Ad jacen t b lock

2 實(shí)驗(yàn)與討論

2.1 試驗(yàn)數(shù)據(jù)及濾波結(jié)果

本文選取國(guó)際攝影測(cè)量與遙感學(xué)會(huì) (ISPRS)在線(xiàn)發(fā)布的參考數(shù)據(jù)作為試驗(yàn)數(shù)據(jù)。對(duì)每個(gè)參考數(shù)據(jù)都事先做好了人工分類(lèi),標(biāo)上地面點(diǎn)或者非地面點(diǎn),以利于檢驗(yàn)濾波的正確性。共選取了 6套測(cè)試數(shù)據(jù),這些數(shù)據(jù)中既有城區(qū)地形也有農(nóng)村地形,其中測(cè)區(qū) 1(Sample 24)的典型地物為陡坎、人字型房屋和植被,測(cè)區(qū) 2(Sample 31)的典型地物為一個(gè)方形復(fù)雜的建筑物,測(cè)區(qū) 3(Sample 12)的典型地物為城區(qū)密集建筑,測(cè)區(qū) 4(Sample 41)的典型地物為形狀不規(guī)則層次多的建筑,測(cè)區(qū) 5(Sample 51)的典型地物為斜坡及斜坡上的植物,測(cè)區(qū) 6(Sample 54)為密集低矮植物。圖 6顯示了測(cè)區(qū) 3和測(cè)區(qū) 6的原始數(shù)據(jù)及其濾波結(jié)果。

圖 6 測(cè)試數(shù)據(jù)與濾波結(jié)果Fig.6 Raw data and filtering results

2.2 濾波結(jié)果分析

從圖 6可以看出,該算法能夠較好地濾除密集的規(guī)則房屋 (測(cè)區(qū) 3)等人工建筑物,對(duì)于平地上低矮的植物和斜坡 (測(cè)區(qū) 6)也有較好的效果。但是,該算法對(duì)于地表斷裂處 (測(cè)區(qū) 1)存在一定的誤差。

表 1列出 6個(gè)測(cè)區(qū)的原方法與改進(jìn)方法的濾波結(jié)果,包括每個(gè)測(cè)區(qū)中標(biāo)準(zhǔn)數(shù)據(jù)和兩種濾波方法濾出的正確地面點(diǎn)及地物點(diǎn)點(diǎn)數(shù)。

表 1 各測(cè)區(qū)濾波結(jié)果比較Tab.1 Comparasion of filtering results (個(gè))

從表 1可以看出,線(xiàn)性迭代最小二乘濾波可以較好地濾出地物點(diǎn),這樣可以保證得到的DEM有較高的精度;但有時(shí)原算法的地面點(diǎn)分類(lèi)正確率相對(duì)較低 (如測(cè)區(qū) 4),而改進(jìn)的方法在盡量保證地物點(diǎn)濾波正確性的同時(shí),提高了地面點(diǎn)的正確率。

圖 7示出 6個(gè)測(cè)區(qū)的誤差統(tǒng)計(jì)圖。假定 a表示數(shù)據(jù)中地面點(diǎn)的點(diǎn)數(shù),b表示數(shù)據(jù)中地物點(diǎn)點(diǎn)數(shù),c表示分類(lèi)正確的地面點(diǎn)點(diǎn)數(shù),d表示分類(lèi)正確的地物點(diǎn)點(diǎn)數(shù),則正確分類(lèi)的地面點(diǎn)比例為 100×%,正確分類(lèi)的地物點(diǎn)的比例為 100×%,總分類(lèi)正確點(diǎn)的比例為 100×%。

圖 7 濾波結(jié)果正確率比較Fig.7 Comparison of filter ing accuracy

從圖 7(a)可以看出,原濾波算法的總點(diǎn)數(shù)精度在 90%左右,但波動(dòng)較大,說(shuō)明該方法對(duì)于不同的地形有著一定的依賴(lài)性。例如對(duì)于含有較密集建筑的第 3測(cè)區(qū),濾波結(jié)果就不如濾除植被的第 6測(cè)區(qū)好。而改進(jìn)的方法,由于提前對(duì)特殊地形進(jìn)行了分析,在處理房屋時(shí)就提高了濾波精度。圖 7(b)改進(jìn)的方法比原方法的地面正確率更高。從單個(gè)測(cè)區(qū)看 ,第 2、4、6測(cè)區(qū)的正確率高于第 1、3、5測(cè)區(qū) ,說(shuō)明該方法在陡坎、斜坡密集的地方,濾出的地面點(diǎn)精度會(huì)降低。圖 7(c)除了第 1測(cè)區(qū)外,其他測(cè)區(qū)精度都較高。由于第 1測(cè)區(qū)中有一個(gè)陡坎,陡坎有高程突變,濾波時(shí)會(huì)誤檢測(cè)為地物點(diǎn),造成一定誤差。

總的來(lái)看,無(wú)論是城區(qū)的復(fù)雜建筑、密集房屋,還是農(nóng)村的斜坡、低矮植被,改進(jìn)的濾波算法都有較好的濾波結(jié)果。從地面點(diǎn)和地物點(diǎn)分類(lèi)正確率來(lái)看,改進(jìn)的算法在盡量保證地物分類(lèi)正確的同時(shí),提高了地面點(diǎn)分類(lèi)的正確率。

3 結(jié)論

(1)在總結(jié)各種濾波算法基礎(chǔ)上,通過(guò)對(duì)高程直方圖的分析,在粗差剔除、初始點(diǎn)選取、特殊地形分析方面對(duì)線(xiàn)性迭代最小二乘濾波方法進(jìn)行了改進(jìn),并在濾波結(jié)束后檢驗(yàn)了相鄰分塊的光滑程度。

(2)采用總點(diǎn)數(shù)分類(lèi)正確率、地面點(diǎn)正確率和地物點(diǎn)正確率 3個(gè)評(píng)價(jià)指標(biāo)對(duì)改進(jìn)的方法和原線(xiàn)性迭代最小二乘濾波結(jié)果進(jìn)行了比較。結(jié)果表明,改進(jìn)的方法在處理密集房屋等方面提高了濾波精度。

[1] Vosselman S G.Experimental Comparison of Filter Algorithms for Bare-Earth Extraction from Airborne Laser Scanning Point Clouds[J].ISPRS Journal of Photogrammetry and Remote Sensing,2004,59(1/2):85-101.

[2] Lindenberger J.Laser-profilmenssungen zur Topographischen Gelandeaufnahm e[D].Stuttgart:Universitat Stuttgart,Verleg der Bayerischen Akademie der Wissenschaften,1993.

[3] Kilian J,Haala N,EnglichM.Cap ture and Evaluation of Airborne Laser ScannerData[J].IAPRS,1996(31-B3):383.

[4] Zhang K Q,Chen S,Whitman D,et al.Progressive Morphological Filter for Removing Nonground Measurements from Airborne L i-DAR Data[J].IEEE Transactions on Geoscience and Remote Sensing,2003,41(4):872-882.

[5] Axelsson P.DEM Generation from Laser Scanner Data Using Adaptive TIN Models[J].International Archives of the Photogrammetry,Remote Sensing and Spatial Information Sciences,2000,33(B4/1):110-117.

[6] Kraus K,Pfeifer N.Determination of Terrain Models in Wooded Areas with Airborne Laser Scanner Data[J].ISPRS Journal of Photogrammetry and Remote Sensing,1998,54(4):193-203.

[7] Kraus K,Pfeifer N.Advance DTM Generation from LiDAR Data[J].International Archives of Photogrammetry and Remote Sensing,2001,34(3/W 4):23-35.

[8] 劉經(jīng)南,許曉東,張小紅,等.機(jī)載激光掃描測(cè)高數(shù)據(jù)分層迭代選權(quán)濾波方法及其質(zhì)量評(píng)價(jià)[J].武漢大學(xué)學(xué)報(bào) (信息科學(xué)版),2008,33(6):551-555.

[9] 王 刃,徐 青,朱新慧,等.用分層穩(wěn)健線(xiàn)性估計(jì)法從機(jī)載 L i-DAR數(shù)據(jù)中獲取DEM[J].測(cè)繪科學(xué)技術(shù)學(xué)報(bào),2008,25(3):188-191.

[10]曾齊紅,毛建華,李先華,等.機(jī)載激光雷達(dá)點(diǎn)云的階層式分類(lèi)[J].測(cè)繪科學(xué),2008,33(1):103-105.

(責(zé)任編輯:劉心季)

Progressive Linear Prediction Fitting for Extracting DTM from Airborne LiDAR Data

ZHANG Jing1,ZHANG Xiao-jun1,JIANGW an-shou1,WANG Jian-chao2,GUO Da-hai2
(1.LIESMARS,Wuhan University,Wuhan 430079,China;2.China Aero Geophysical Survey&Remote Sensing Center for Land and Resousces,Beijing 100083,China)

A progressive linear prediction filtering algorithm is proposed for extracting DEM from LiDAR data.Some processes are inserted in the ordinary linear prediction algorithm,such as gross error detection,initial value selection,land form analysis,and smoothness detection.The authors used this algorithm to process 6 datasets published by ISPRS as standard filtering test data.The results show that the improvement in the traditional methods can increase the precision of DEM.

LiDAR;Filtering;Linear prediction;Outlier;Partitioning

江萬(wàn)壽 (1967-),男,研究員,博士生導(dǎo)師。主要研究方向?yàn)閿?shù)字?jǐn)z影測(cè)量及遙感數(shù)據(jù)處理。聯(lián)系電話(huà):(+86)27-68778092-8321(O),E-mail:jw s@lm ars.whu.edu.cn。

TP 75

A

1001-070X(2011)01-0052-05

2010-04-21;

2010-06-10

國(guó)家“863”計(jì)劃項(xiàng)目 (編號(hào):2006AA 06A 208)和國(guó)家自然科學(xué)基金項(xiàng)目 (編號(hào):40671159)共同資助。

張 靖 (1982-),男,博士研究生,主要從事 LiDAR數(shù)據(jù)處理方面的研究。