楊宇妍,臧玉府,肖雄武,管海燕,彭代峰

1. 南京信息工程大學(xué)遙感與測(cè)繪工程學(xué)院,江蘇 南京 210044; 2. 武漢大學(xué)測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,湖北 武漢 430079

機(jī)載激光掃描系統(tǒng)可高效、準(zhǔn)確地獲取大范圍區(qū)域的空間點(diǎn)云數(shù)據(jù),其中地面點(diǎn)是建立高精度數(shù)字高程模型的重要數(shù)據(jù)源,已在地形地貌、水土保持、土壤侵蝕分析等領(lǐng)域中得到廣泛應(yīng)用[1]。機(jī)載激光點(diǎn)云濾波作為識(shí)別地面點(diǎn)的主要手段,對(duì)地面信息提取、道路規(guī)劃及目標(biāo)精準(zhǔn)識(shí)別等應(yīng)用具有重要作用。但由于點(diǎn)云數(shù)量巨大,且地物易被遮擋,使得處理效率低下,目標(biāo)分類(lèi)困難,對(duì)現(xiàn)有濾波方法提出了挑戰(zhàn)[2]。

目前,機(jī)載激光點(diǎn)云的濾波方法主要包括4類(lèi)[3]:①基于坡度變化的方法。逐點(diǎn)計(jì)算點(diǎn)與相鄰點(diǎn)之間的坡度,通過(guò)比較兩點(diǎn)間實(shí)際坡度與坡度閾值濾除地物點(diǎn)[4-6]。但這類(lèi)方法在處理陡坡、高程突變區(qū)域時(shí)會(huì)錯(cuò)誤濾除大量地面點(diǎn),且坡度越陡濾波精度越低。②數(shù)學(xué)形態(tài)學(xué)方法。對(duì)地形進(jìn)行腐蝕膨脹,構(gòu)造開(kāi)運(yùn)算,利用開(kāi)運(yùn)算后非地面點(diǎn)高程變化大這一特點(diǎn),將高差大于閾值的點(diǎn)歸為非地面點(diǎn)并進(jìn)行濾除[7-9]。這類(lèi)方法得到的濾波結(jié)果依賴(lài)于最佳窗口尺寸的選擇,較小的窗戶(hù)尺寸會(huì)保留較大的建筑物,較大的尺寸會(huì)忽視地形細(xì)節(jié)。③基于曲面的方法。通過(guò)迭代從原始數(shù)據(jù)集中選擇地面測(cè)量值,逐漸逼近地面,創(chuàng)建一個(gè)近似裸露的地球表面[1,9-13]。文獻(xiàn)[12]利用種子點(diǎn)創(chuàng)建稀疏三角網(wǎng),并在迭代過(guò)程中加密,用來(lái)處理不連續(xù)的表面。文獻(xiàn)[13]提出基于漸進(jìn)三角網(wǎng)濾波方法,它在經(jīng)典三角網(wǎng)濾波的基礎(chǔ)上,獲取更具可靠性的種子點(diǎn),以增強(qiáng)對(duì)不同地形的適應(yīng)性。不同于傳統(tǒng)基于曲面的方法,布料模擬濾波(cloth simulation filter,CSF)[14]將點(diǎn)云翻轉(zhuǎn),假設(shè)有一塊布料受到重力落下,最終落下的布代表當(dāng)前地形。此類(lèi)方法在地形變化平緩的地區(qū)效果較好、所需設(shè)置的參數(shù)少。④基于采樣線(xiàn)的濾波方法。此類(lèi)方法因能顧及整個(gè)區(qū)域的地形特征,高效地處理大范圍場(chǎng)景區(qū)域,為機(jī)載激光點(diǎn)云濾波提供了一種思路。文獻(xiàn)[15]中用球面坐標(biāo)代替笛卡兒坐標(biāo),利用方位角和徑向距離球面坐標(biāo)作為判斷標(biāo)準(zhǔn),通過(guò)球坐標(biāo)中徑向距離曲線(xiàn)的斷點(diǎn)和轉(zhuǎn)折點(diǎn)分割地面點(diǎn)和非地面點(diǎn),可有效避免過(guò)度分割和欠分割問(wèn)題。文獻(xiàn)[16]針對(duì)每條采樣線(xiàn),利用半徑可變的圓環(huán)滾過(guò),保留采樣線(xiàn)上被圓環(huán)滾過(guò)的點(diǎn),采用B樣條擬合曲面,通過(guò)比較擬合高程與實(shí)際高程的差值得到地面點(diǎn)。文獻(xiàn)[17]提出一種基于采樣線(xiàn)處理的地面濾波算法,利用一維迭代樣條插值算法改進(jìn)地面點(diǎn)的提取。

現(xiàn)有濾波模型雖能基于一維地形特征實(shí)現(xiàn)快速濾波,但不能避免斷裂線(xiàn)周?chē)拇罅康孛纥c(diǎn)被錯(cuò)誤濾除,難以適用于由斷裂線(xiàn)引起的高程突變區(qū)。針對(duì)以上問(wèn)題,本文基于采樣線(xiàn)濾波優(yōu)勢(shì),構(gòu)建融合斷裂線(xiàn)約束的采樣線(xiàn)濾波模型,使其在濾波過(guò)程中較大程度保留地形斷裂處的地面點(diǎn),適用于多種復(fù)雜地面場(chǎng)景的地物點(diǎn)濾除,提高濾波的完整度和準(zhǔn)確度。

1 基于斷裂線(xiàn)的激光點(diǎn)云濾波

本文通過(guò)建立斷裂線(xiàn)約束模型提出一種高效的采樣線(xiàn)濾波方法。首先,基于點(diǎn)曲率及鄰域范圍內(nèi)距離之和創(chuàng)建約束能量項(xiàng),改進(jìn)Snake能量模型[18],并通過(guò)能量函數(shù)最小化提取斷裂線(xiàn)特征;然后,通過(guò)構(gòu)建正三角形模型,計(jì)算點(diǎn)云平坦度及融合斷裂線(xiàn)約束進(jìn)行采樣線(xiàn)濾波;最后,根據(jù)采樣線(xiàn)上的地面點(diǎn)利用改進(jìn)最小二乘法擬合全局曲面,精準(zhǔn)濾除地物點(diǎn)。該方法主要流程如圖1所示。

1.1 地形斷裂處點(diǎn)云濾波問(wèn)題

采集區(qū)域內(nèi)通常場(chǎng)景復(fù)雜,存在大量的高程突變區(qū)域,嚴(yán)重影響了現(xiàn)有濾波方法的性能,一直是點(diǎn)云濾波中的重難點(diǎn)。研究發(fā)現(xiàn)點(diǎn)云中高程突變區(qū)域附近往往存在顯著的斷裂線(xiàn)特征,為充分利用該信息,本文按地形表面是否連續(xù)的特點(diǎn),將斷裂線(xiàn)特征分為以下兩類(lèi):山谷、山脊等處地形表面連續(xù)變化且存在明顯線(xiàn)特征的Ⅰ類(lèi)斷裂線(xiàn)(圖2(a)),以及在人工建筑、陡坎、沖溝等地表不連續(xù)且存在明顯高程突變的Ⅱ類(lèi)斷裂線(xiàn)(圖2(b))。圖2中黃色曲線(xiàn)為地形斷裂線(xiàn)特征?，F(xiàn)有濾波方法大多沒(méi)有綜合考慮上述地形斷裂線(xiàn)的影響,使得地形突變區(qū)域的地面點(diǎn)或地物點(diǎn)被錯(cuò)誤地濾除或保留。針對(duì)上述問(wèn)題,本文提出一種融合斷裂線(xiàn)約束模型的采樣線(xiàn)高效濾波方法,提高濾波精度、避免高程突變區(qū)域的地形失真。

圖2 斷裂線(xiàn)Fig.2 Breaklines

1.2 地形斷裂線(xiàn)提取

曲率是反映地形表面起伏程度的重要度量,斷裂線(xiàn)上的點(diǎn)大多具有局部曲率極值。本文基于點(diǎn)曲率及鄰域范圍內(nèi)距離之和創(chuàng)建約束能量項(xiàng),改進(jìn)Snake能量模型,并通過(guò)能量函數(shù)最小化提取斷裂線(xiàn)特征。對(duì)點(diǎn)云中任意一點(diǎn)及其鄰域點(diǎn),采用最小二乘擬合法實(shí)現(xiàn)其局部曲面參數(shù)化,設(shè)該點(diǎn)局部曲面函數(shù)為f(xsur,ysur),其一階、二階偏導(dǎo)數(shù)分別表示為rx、ry、rxx、rxy、ryy,則該點(diǎn)處的平均曲率H(s)[19]為

(1)

(2)

式中,Eint表示模型內(nèi)部能量,控制曲線(xiàn)的平滑性和連續(xù)性;Eext表示模型外部能量,吸引初始線(xiàn)向著曲線(xiàn)延伸或伸縮;Econ表示模型約束能量,使初始斷裂線(xiàn)向外部能量較弱的地方繼續(xù)演化。內(nèi)部能量Eint、外部能量Eext和約束能量Econ構(gòu)建如下

(3)

Eext[l(s)]=ω(s)·|u(s)|

(4)

(5)

通過(guò)循環(huán)迭代求解能量函數(shù)最小值[22],排除噪聲點(diǎn)的影響、優(yōu)化Snake曲線(xiàn),得到最終的斷裂線(xiàn),如圖3所示。

圖3 部分?jǐn)嗔丫€(xiàn)提取結(jié)果Fig.3 Results of partial breaklines extraction

1.3 融合斷裂線(xiàn)約束的濾波模型

為滿(mǎn)足大范圍點(diǎn)云濾波需求,本文由測(cè)區(qū)的機(jī)載點(diǎn)云生成多條采樣線(xiàn),將空間濾波問(wèn)題轉(zhuǎn)換為一維濾波處理。在點(diǎn)云場(chǎng)景范圍內(nèi)沿縱坐標(biāo)方向等間距對(duì)點(diǎn)云進(jìn)行采樣,獲取點(diǎn)序列構(gòu)成的采樣線(xiàn)(圖4),再將每根采樣線(xiàn)沿縱軸方向上給定距離d內(nèi)的點(diǎn)歸入對(duì)應(yīng)的采樣線(xiàn),將采樣線(xiàn)上的點(diǎn)按橫坐標(biāo)值升序排列。距離d的大小會(huì)影響采樣線(xiàn)準(zhǔn)確性:若給定距離過(guò)大,則采樣線(xiàn)上的點(diǎn)數(shù)量增多,造成數(shù)據(jù)冗余,繼而影響采樣線(xiàn)濾波的效率;反之,采樣線(xiàn)上的點(diǎn)數(shù)量過(guò)少,采樣線(xiàn)喪失連續(xù)性,不但降低采樣線(xiàn)濾波精度,同時(shí)也使得曲面擬合喪失真實(shí)性,故本文通過(guò)式(6)計(jì)算距離d[15]

(6)

圖4 采樣線(xiàn)Fig.4 Sampling line

式中,num為點(diǎn)云數(shù)量;S和W分別表示點(diǎn)云所表示地形范圍的長(zhǎng)和寬。

由圖4可知,本文求得的采樣線(xiàn)連貫、準(zhǔn)確,能保證后續(xù)操作順利進(jìn)行。

1.3.1 預(yù)處理

通過(guò)構(gòu)建等邊三角形模型、計(jì)算點(diǎn)平坦度及鄰域范圍內(nèi)高差進(jìn)行采樣線(xiàn)粗濾波(圖5),濾除采樣線(xiàn)上的車(chē)輛、植被、建筑物等非地面點(diǎn),再選取粗濾波后采樣線(xiàn)上的最左端點(diǎn)為采樣中心點(diǎn),取距采樣中心點(diǎn)最近的點(diǎn)為基準(zhǔn)點(diǎn)。

圖5 采樣線(xiàn)粗過(guò)濾Fig.5 Coarse filtering of a sampling line

(7)

(8)

在△O1P2O2中,運(yùn)用三角余弦定理計(jì)算O1P2和P2O2的夾角θ,判斷θ與夾角閾值θthreshold的關(guān)系:若θ大于θthreshold,則P3為非地面點(diǎn),繼續(xù)取采樣線(xiàn)上與P2相鄰下一點(diǎn)作為目標(biāo)點(diǎn),重復(fù)上述步驟;反之,P3為地面點(diǎn),用P2更新P1、P3更新P2繼續(xù)判斷后續(xù)點(diǎn),直至采樣線(xiàn)終點(diǎn)。

濾除采樣線(xiàn)上大部分非地面點(diǎn)后,需要利用點(diǎn)云平坦度λ[23]進(jìn)一步過(guò)濾植被、車(chē)輛和人員等地物。本文根據(jù)上述處理后采樣線(xiàn)點(diǎn)局部鄰域R內(nèi)原始點(diǎn)Xi(0

(9)

(10)

計(jì)算點(diǎn)云平坦度濾除采樣線(xiàn)上的植被、車(chē)輛和人員等地物后,需要考慮建筑物等局部鄰域平坦的非地面點(diǎn)對(duì)采樣中心選取的影響[24]。為剔除采樣線(xiàn)上的建筑物點(diǎn),對(duì)點(diǎn)云中任一點(diǎn),計(jì)算統(tǒng)計(jì)高程大于該點(diǎn)的鄰域點(diǎn)數(shù)與全部鄰域點(diǎn)數(shù)的比值,比值越小,該點(diǎn)為地面點(diǎn)的可能性越大。對(duì)每條采樣線(xiàn)進(jìn)行粗過(guò)濾后,采樣線(xiàn)上最左端點(diǎn)為采樣中心,距采樣中心最近的點(diǎn)為基準(zhǔn)點(diǎn)。

1.3.2 基于斷裂線(xiàn)的采樣線(xiàn)點(diǎn)云遍歷

對(duì)任一條采樣線(xiàn)精細(xì)濾波,基于斷裂線(xiàn)約束,結(jié)合采樣中心點(diǎn)、基準(zhǔn)點(diǎn)和各采樣點(diǎn)構(gòu)建有向向量計(jì)算夾角[25],濾除地物點(diǎn)。確定采樣中心點(diǎn)和基準(zhǔn)點(diǎn)后,從原始采樣線(xiàn)最左端開(kāi)始遍歷,由近及遠(yuǎn)按點(diǎn)云順序計(jì)算坡度ρ,得到采樣線(xiàn)上的地面點(diǎn)。坡度ρ的計(jì)算公式為

(11)

圖6為基于斷裂線(xiàn)的采樣線(xiàn)點(diǎn)云遍歷示意圖,具體的遍歷步驟如下。

圖6 采樣線(xiàn)點(diǎn)云遍歷Fig.6 A sampling line point cloud traversal

∥輸入:采樣中心點(diǎn)O、基準(zhǔn)點(diǎn)B、采樣點(diǎn)數(shù)n

fori←0 ton

ifPi==斷裂點(diǎn)

保留Pi

ifPi!=斷裂點(diǎn)

連接BPi、OB,根據(jù)式(11)計(jì)算OB和BPi的夾角ρ

ifρ≤ρthreshold(閾值)

保留Pi,Pi=B,i++

else

濾除Pi,i++

∥輸出:一維地面點(diǎn)(GRPts)

基于斷裂線(xiàn)約束的采樣線(xiàn)濾波結(jié)果如圖7所示,由圖7(b)可知,本文方法能準(zhǔn)確區(qū)分采樣線(xiàn)上藍(lán)色地面點(diǎn)與紅色非地面點(diǎn),確保后續(xù)處理順利執(zhí)行。

圖7 采樣線(xiàn)濾波前后對(duì)比結(jié)果Fig.7 Comparison before and after sampling lines filtering

1.4 改進(jìn)的二次曲面擬合

由于地表復(fù)雜多樣,簡(jiǎn)單的二次曲面擬合方法不能真實(shí)表示濾波后的地表形態(tài),在連續(xù)性上存在較大缺陷。為較大程度還原地面,保證地表的真實(shí)性和連續(xù)性,本文格網(wǎng)化整個(gè)測(cè)區(qū)(格網(wǎng)單元尺寸原則上應(yīng)大于最大地物的邊長(zhǎng)),且各格網(wǎng)的邊界區(qū)域相交重疊(圖8)。圖8中黃色區(qū)域表示任一格網(wǎng)的不重疊區(qū)域,綠色區(qū)域表示兩個(gè)格網(wǎng)的重疊部分,藍(lán)色區(qū)域表示4個(gè)格網(wǎng)的重疊部分。先根據(jù)落在單元格(圖8中虛線(xiàn)矩形)內(nèi)的處理后的m個(gè)采樣線(xiàn)上的地面點(diǎn),運(yùn)用最小二乘原理計(jì)算誤差平方和σ2達(dá)到最小時(shí)的近似地面參數(shù)A、B、C、D、E、F,依據(jù)地面參數(shù)計(jì)算點(diǎn)云的擬合高程z(x,y)。σ2和z(x,y)的計(jì)算公式為

圖8 格網(wǎng)劃分Fig.8 Grids division

(12)

z(x,y)=Ax2+By2+Cxy+Dx+Ey+F

(13)

式中,σ2為計(jì)算值與真實(shí)值之間的誤差平方和;z(x,y)為點(diǎn)云在二次曲面上的高程。判斷原始點(diǎn)云中每個(gè)點(diǎn)所在的格網(wǎng)位置:若原始點(diǎn)落入非重疊區(qū)域,則依據(jù)式(13)計(jì)算點(diǎn)云的擬合高程z(x,y);若點(diǎn)落入格網(wǎng)重疊區(qū)域,則計(jì)算點(diǎn)在不同格網(wǎng)曲面擬合高程z(x,y)的平均值。比較擬合值與實(shí)際高程的差值,將差值大于閾值的點(diǎn)從候選地形點(diǎn)中去除。

2 試驗(yàn)結(jié)果與分析

為驗(yàn)證本文方法的有效性,采用國(guó)際攝影測(cè)量與遙感學(xué)會(huì)(ISPRS)工作組提供的數(shù)據(jù)集和國(guó)內(nèi)成都山區(qū)點(diǎn)云兩套數(shù)據(jù)進(jìn)行測(cè)試。參考數(shù)據(jù)集是通過(guò)手動(dòng)過(guò)濾點(diǎn)云數(shù)據(jù)集生成的,地形樣本中的每個(gè)點(diǎn)被分類(lèi)為地面點(diǎn)或地物點(diǎn),后續(xù)將手動(dòng)分類(lèi)結(jié)果作為真值檢驗(yàn)方法精度。

2.1 ISPRS試驗(yàn)數(shù)據(jù)

對(duì)于ISPRS提供的數(shù)據(jù),本文選擇了4組具有不同地形特征的數(shù)據(jù)集對(duì)方法性能進(jìn)行測(cè)試,詳細(xì)信息見(jiàn)表1。

表1 地形特征

2.1.1 試驗(yàn)結(jié)果

為了評(píng)估方法的性能,本文計(jì)算了4個(gè)樣本的Ⅰ型誤差、Ⅱ型誤差和總誤差[1]。此外,還通過(guò)Kappa系數(shù)[26]衡量分類(lèi)結(jié)果的總體一致性。各型誤差、Kappa系數(shù)和處理時(shí)間見(jiàn)表2。結(jié)果表明,對(duì)于不同的地形特征,本文方法均能有效提取地面點(diǎn),并將總濾波誤差控制在5%左右。其中地形2的總誤差最小,這表明本文方法在存在大型建筑和橋梁、且地形平坦的情形下,性能依然優(yōu)異。本文方法在處理Ⅰ型誤差時(shí),表現(xiàn)尤為突出,均控制在2%左右,這表明該方法能較大程度地保留地面點(diǎn),還原地面形態(tài)。4個(gè)地形的濾波時(shí)間分別為22.6、18.3、25.4、8.3 s,可見(jiàn)本文方法兼顧了精度和效率。總體而言,本文方法對(duì)地形的起伏較為敏感,能夠處理各種地貌類(lèi)型,并且在斜坡、建筑邊緣起伏處表現(xiàn)良好。

表2 場(chǎng)景誤差分析

試驗(yàn)中4個(gè)典型地形濾波結(jié)果如圖9所示。圖9(a)、(c)、(e)、(g)為濾波前地形,圖9(b)、(d)、(f)、(h)展示了各類(lèi)地形中附著在DEM表面上的Ⅱ型誤差(紅色點(diǎn)),DEM由地面點(diǎn)生成。從圖9中對(duì)不同地形的濾波結(jié)果看,疏密程度不同的地表植被和各種形狀的建筑物基本被濾除,同時(shí)保留了地形特征和地面細(xì)節(jié)。

圖9 濾波結(jié)果Fig.9 Filtering results

為了全面分析本文方法的準(zhǔn)確性,將本文方法的總誤差和Kappa系數(shù)與現(xiàn)有方法[9,14,27-29]進(jìn)行了比較。這些方法和本文方法的總誤差與Kappa系數(shù)如圖10—圖11所示?？傮w而言,對(duì)于不同特征的地形,本文方法通用性好并且濾波誤差更小。除地形3外,本文方法與傳統(tǒng)方法相比濾波精度提高了3%～5%。Kappa系數(shù)是統(tǒng)計(jì)學(xué)中度量一致性的指標(biāo)。Kappa系數(shù)越大,表示模型預(yù)測(cè)結(jié)果與實(shí)際濾波結(jié)果一致性越強(qiáng);若Kappa系數(shù)在[0.81,1]區(qū)間,則說(shuō)明模擬結(jié)果與實(shí)際結(jié)果幾乎完全一致,本文方法計(jì)算得到的Kappa系數(shù)均高于88%?？傉`差和Kappa系數(shù)的比較結(jié)果,說(shuō)明本文方法準(zhǔn)確性、精度與現(xiàn)有方法比較均表現(xiàn)優(yōu)異。

圖10 本文方法與其他方法的總誤差比較Fig.10 Total error compared between the proposed method and other reported methods

圖11 本文方法與其他方法的Kappa系數(shù)比較Fig.11 Kappa coefficient compared between the proposed method and other reported methods

2.1.2 參數(shù)設(shè)置

本文將通過(guò)大量測(cè)試確定的普適參數(shù)設(shè)置為固定值,在此基礎(chǔ)上,針對(duì)不同樣本確定采樣線(xiàn)間距、坡度、格網(wǎng)邊長(zhǎng)及曲率計(jì)算時(shí)的鄰域范圍等參數(shù)的最優(yōu)值。為了定量評(píng)估采樣線(xiàn)間距的影響,本文測(cè)試了具有不同間距的樣本(從6～20 m,步長(zhǎng)為2 m)。每組總誤差(取決采樣線(xiàn)間距)如圖12所示。數(shù)據(jù)表明,所有組的總誤差基本隨著采樣線(xiàn)間距增大而增大,只有地形1在間距最大(20 m)的時(shí)候,取得最好的濾波效果。若采樣線(xiàn)間距太小,會(huì)得到大量采樣線(xiàn)點(diǎn)云,將失去簡(jiǎn)化數(shù)據(jù)、節(jié)省時(shí)間的意義。若采樣線(xiàn)間距太大,后續(xù)生成的格網(wǎng)中可能不存在處理后的地面點(diǎn),導(dǎo)致曲面擬合失敗,無(wú)法對(duì)全域進(jìn)行濾波。因此,大多樣本均選擇6 m作為采樣線(xiàn)間距值,因?yàn)樗鼘?duì)所有樣本都產(chǎn)生了相對(duì)較好的結(jié)果,并且可以應(yīng)用于許多情況而無(wú)須調(diào)整。

濾波過(guò)程中采樣線(xiàn)上地物點(diǎn)濾除效果與坡度值設(shè)置有較強(qiáng)的關(guān)系。圖13顯示了不同坡度值下的總誤差變化,可以看出,隨著坡度的增加,濾波總誤差也隨之增大。所以所有樣本的最佳坡度值均為20°。

圖13 坡度閾值對(duì)應(yīng)的總誤差Fig.13 Total errors for each slope threshold

曲面擬合過(guò)程中格網(wǎng)邊長(zhǎng)將影響擬合效果(圖14)。如圖14所示,若格網(wǎng)邊長(zhǎng)設(shè)置過(guò)小,格網(wǎng)內(nèi)采樣線(xiàn)上的地面點(diǎn)數(shù)可能為0,將無(wú)法求解多項(xiàng)式系數(shù),導(dǎo)致擬合曲面失敗;若格網(wǎng)邊長(zhǎng)設(shè)置過(guò)大,擬合的曲面將不能反映真實(shí)的地表情況。因此,大多地形選擇邊長(zhǎng)為40 m的格網(wǎng),此時(shí)的擬合效果最佳。

圖14 格網(wǎng)邊長(zhǎng)對(duì)應(yīng)的總誤差Fig.14 Total errors for each grid size

除上述參數(shù)外,斷裂線(xiàn)的完整性與濾波結(jié)果緊密相關(guān)。本文采用改進(jìn)Snake模型提取斷裂線(xiàn),曲率作用于約束能量,故計(jì)算曲率時(shí)鄰域的選取制約著斷裂線(xiàn)的提取(圖15)。由圖15可知,若鄰域范圍選取過(guò)大,造成斷裂線(xiàn)冗余,將會(huì)影響采樣線(xiàn)濾波結(jié)果,從而導(dǎo)致擬合曲面過(guò)高、錯(cuò)誤濾除地面點(diǎn);相反,斷裂線(xiàn)提取不完整將不能濾除高程突變區(qū)的地物點(diǎn)。

圖15 曲率計(jì)算時(shí)鄰域范圍對(duì)應(yīng)的總誤差Fig.15 Total error corresponding to different neighborhood in curvature calculation

2.1.3 消融研究

對(duì)于場(chǎng)景中包含許多陡坡或建筑的地形,斷裂線(xiàn)是影響本文方法準(zhǔn)確性的重要因素。假設(shè)不考慮斷裂線(xiàn)對(duì)濾波結(jié)果的影響,那么在陡坡、建筑物邊緣處的地面點(diǎn)將被錯(cuò)誤地移除,影響最終的濾波精度。解決此問(wèn)題的一種方法是提取整塊區(qū)域的斷裂線(xiàn),并在隨后的采樣線(xiàn)濾波中保留斷裂點(diǎn),正確處理陡坡區(qū)域。加入斷裂線(xiàn)約束前后濾波結(jié)果見(jiàn)表3,可以看出4類(lèi)地形的總誤差均減小,其中地形2的總誤差減少近10%;Ⅰ型誤差也大幅減小,這表明在加入斷裂線(xiàn)約束后,在陡峭地表邊緣的地面點(diǎn)會(huì)被正確保留;與此同時(shí)緊貼地表的地物點(diǎn)會(huì)被視為地面點(diǎn)而被錯(cuò)誤保留,因此Ⅱ型誤差略有增大。

表3 添加斷裂線(xiàn)約束前后的誤差比較

傳統(tǒng)Snake模型利用內(nèi)部能量控制曲線(xiàn)的平滑性和連續(xù)性,利用外部能量限制Snake的演變,而初始輪廓線(xiàn)距斷裂線(xiàn)較遠(yuǎn)時(shí),外部能量的作用不明顯,故本文加入曲率約束和距離約束構(gòu)成能量約束項(xiàng),用于輔助外部能量對(duì)Snake演變的限制。Snake模型改進(jìn)前后的濾波結(jié)果見(jiàn)表4,可以看出,采用傳統(tǒng)的Snake模型提取斷裂線(xiàn)能保證濾波精度在90%以上,這表明本文基于斷裂線(xiàn)約束的濾波方法頗具成效。在Snake模型中加入約束能量項(xiàng)后,地形2的濾波精度能達(dá)到97.18%,其余地形的總誤差也控制在5%左右,表明改進(jìn)Snake模型的有效性。

采用傳統(tǒng)二次曲面擬合方法對(duì)測(cè)區(qū)進(jìn)行整體擬合,不能較大程度逼近真實(shí)地面,故本文按測(cè)區(qū)點(diǎn)云的x、y坐標(biāo)劃分格網(wǎng),進(jìn)行格網(wǎng)局部擬合,且各格網(wǎng)的邊界區(qū)域相交重疊,提高了濾波精度。二次曲面擬合方法改進(jìn)前后的濾波結(jié)果見(jiàn)表5。由表5可知,除地形3的Ⅱ型誤差在改進(jìn)二次曲面擬合后略微增加2%外,其余地形的各類(lèi)誤差均不同程度降低,其中地形4的濾波效果改善最為明顯,濾波精度提高了8.64%,表明改進(jìn)二次曲面擬合方法能真實(shí)展現(xiàn)地面形態(tài),保證了地表的真實(shí)性和連續(xù)性。

表5 二次曲面擬合改進(jìn)前后的誤差比較

2.2 成都山區(qū)試驗(yàn)數(shù)據(jù)

數(shù)據(jù)來(lái)源于成都市山區(qū),區(qū)域面積約28.9萬(wàn)m2,共約157.5萬(wàn)點(diǎn),海拔為408.55～525.17 m,包括陡坡、茂密植被、梯田、河流和低矮建筑,如圖16所示。圖17顯示了本文方法的濾波結(jié)果,其中掃描線(xiàn)間距為12 m、坡度閾值為20°、格網(wǎng)邊長(zhǎng)為40 m、計(jì)算曲率時(shí)的鄰域范圍為6 m。手動(dòng)檢查表明,本文方法在復(fù)雜區(qū)域表現(xiàn)良好,只有少數(shù)錯(cuò)誤點(diǎn)。由于分割不足和樹(shù)木遮擋,地面點(diǎn)被錯(cuò)分為植被點(diǎn)從而造成少量Ⅰ型誤差,如圖17的矩形A、B、C所示。此外,一些附著地面的低矮植被點(diǎn)被錯(cuò)誤提取,如圖17的矩形D、E所示;加入斷裂線(xiàn)約束后,斷裂線(xiàn)處緊貼地表的地物點(diǎn)被視為地面點(diǎn)而錯(cuò)誤保留,因此Ⅱ型誤差略有增大,如圖17的矩形F所示。

圖16 成都山區(qū)機(jī)載點(diǎn)云數(shù)據(jù)Fig.16 Point clouds from the city of Chengdu

圖17 成都山區(qū)點(diǎn)云濾波結(jié)果Fig.17 The filtering results of the Chengdu data

比較本文方法與商業(yè)軟件Terrasolid TerraScan中布料濾波方法,CC(CloudCompare)中布料濾波、基于坡度濾波方法的各類(lèi)誤差和Kappa系數(shù)(表6),可知本文方法的濾波性能更優(yōu)、表現(xiàn)更好。

表6 本文方法、商業(yè)軟件TerraScan和CloudCompare的誤差比較

圖18顯示了研究區(qū)域的斷裂線(xiàn)及基于斷裂線(xiàn)約束生成的DEM。由圖18中成都山區(qū)的DEM可知,低矮建筑、茂密植被等被濾除,山區(qū)地面的細(xì)節(jié)特征也得以保留。為探究斷裂線(xiàn)對(duì)DEM精度的影響,選取圖18中的3個(gè)區(qū)域(A1、A2、A3)分析有、無(wú)斷裂線(xiàn)約束的DEM與參考地面之間的差異。

圖18 斷裂線(xiàn)提取與DEM生成Fig.18 Extracting brealines and generating DEMs with breaklines

圖19比較了圖18中A1、A2、A3區(qū)域有無(wú)斷裂線(xiàn)約束DEM與參考地面點(diǎn)的二維橫截面,可以看出,具有斷裂線(xiàn)的DEM橫截面更接近參考地面點(diǎn)。此外,由圖19(a)、圖19(b)可知,在無(wú)斷裂線(xiàn)約束時(shí),地形起伏大的地面會(huì)被錯(cuò)誤濾除, 說(shuō)明本文方法在Ⅰ類(lèi)斷裂線(xiàn)和Ⅱ類(lèi)斷裂線(xiàn)處都能準(zhǔn)確地還原地表。在圖19(c)這類(lèi)地形起伏較小的Ⅱ類(lèi)斷裂線(xiàn)處,有無(wú)斷裂線(xiàn)約束都能生成高質(zhì)量DEM,而基于斷裂線(xiàn)約束生成的DEM橫截面更貼合參考地面。

圖19 3個(gè)區(qū)域的有、無(wú)斷裂線(xiàn)約束的DEM二維橫截面比較Fig.19 2D cross-sections comparisons between the DEMs with/without breaklines and the ground truths in three areas

為進(jìn)一步探討斷裂線(xiàn)約束的重要性,本文分別計(jì)算圖18中A1、A2、A3區(qū)域參考DEM與有無(wú)斷裂線(xiàn)約束DEM之間的高差,定量對(duì)比有無(wú)斷裂線(xiàn)約束下的濾波后DEM精度,結(jié)果如圖20所示。誤差主要分布在斷裂線(xiàn)周?chē)?且對(duì)于Ⅰ類(lèi)斷裂線(xiàn)和Ⅱ類(lèi)斷裂線(xiàn)所在區(qū)域,無(wú)斷裂線(xiàn)約束的DEM與參考DEM的高差更大,表明基于斷裂線(xiàn)約束生成的DEM質(zhì)量更高,在地形陡峭處表現(xiàn)更好。

3 結(jié)論與展望

基于地形表面崎嶇不平、存在高程突變的特點(diǎn),本文提出了一種基于斷裂線(xiàn)約束和改進(jìn)二次曲面擬合的點(diǎn)云濾波方法。不同于將單一濾波問(wèn)題轉(zhuǎn)換為一維采樣線(xiàn)問(wèn)題,本文改進(jìn)Snake能量模型,并通過(guò)能量函數(shù)最小化提取斷裂線(xiàn)特征,在采樣線(xiàn)濾波基礎(chǔ)上融合斷裂線(xiàn)約束,阻止了高程突變處的大量地面點(diǎn)被錯(cuò)誤濾除;根據(jù)采樣線(xiàn)上的地面點(diǎn)利用改進(jìn)二次曲面擬合過(guò)濾非地面點(diǎn),較大程度還原地貌,提高了點(diǎn)云濾波精度,增加了點(diǎn)云濾波模型適用的場(chǎng)景類(lèi)型。在4份ISPRS標(biāo)準(zhǔn)數(shù)據(jù)中,濾波的準(zhǔn)確率分別達(dá)到94.74%、97.18%、94.96%和95.26%;在成都山區(qū)數(shù)據(jù)中,濾波精度達(dá)到97.61%。試驗(yàn)結(jié)果表明,本文方法能從機(jī)載激光點(diǎn)云中精確提取地面點(diǎn),不僅適用于平坦的地形表面,而且對(duì)復(fù)雜崎嶇地面也能取得較好的結(jié)果。目前本文濾波模型的參數(shù)設(shè)置還需要一定的先驗(yàn)知識(shí),因此未來(lái)的一個(gè)工作方向是改進(jìn)模型自適應(yīng)參數(shù)選擇,進(jìn)一步提高模型自動(dòng)化程度。