基于地形斷裂線約束的機載激光點云高精度濾波方法

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

1. 南京信息工程大學遙感與測繪工程學院,江蘇 南京 210044; 2. 武漢大學測繪遙感信息工程國家重點實驗室,湖北 武漢 430079

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

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

現有濾波模型雖能基于一維地形特征實現快速濾波,但不能避免斷裂線周圍的大量地面點被錯誤濾除,難以適用于由斷裂線引起的高程突變區。針對以上問題,本文基于采樣線濾波優勢,構建融合斷裂線約束的采樣線濾波模型,使其在濾波過程中較大程度保留地形斷裂處的地面點,適用于多種復雜地面場景的地物點濾除,提高濾波的完整度和準確度。

1 基于斷裂線的激光點云濾波

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

1.1 地形斷裂處點云濾波問題

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

圖2 斷裂線Fig.2 Breaklines

1.2 地形斷裂線提取

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

(1)

(2)

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

(3)

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

(4)

(5)

通過循環迭代求解能量函數最小值[22],排除噪聲點的影響、優化Snake曲線,得到最終的斷裂線,如圖3所示。

圖3 部分斷裂線提取結果Fig.3 Results of partial breaklines extraction

1.3 融合斷裂線約束的濾波模型

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

(6)

圖4 采樣線Fig.4 Sampling line

式中,num為點云數量;S和W分別表示點云所表示地形范圍的長和寬。

由圖4可知,本文求得的采樣線連貫、準確,能保證后續操作順利進行。

1.3.1 預處理

通過構建等邊三角形模型、計算點平坦度及鄰域范圍內高差進行采樣線粗濾波(圖5),濾除采樣線上的車輛、植被、建筑物等非地面點,再選取粗濾波后采樣線上的最左端點為采樣中心點,取距采樣中心點最近的點為基準點。

圖5 采樣線粗過濾Fig.5 Coarse filtering of a sampling line

(7)

(8)

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

濾除采樣線上大部分非地面點后,需要利用點云平坦度λ[23]進一步過濾植被、車輛和人員等地物。本文根據上述處理后采樣線點局部鄰域R內原始點Xi(0

(9)

(10)

計算點云平坦度濾除采樣線上的植被、車輛和人員等地物后,需要考慮建筑物等局部鄰域平坦的非地面點對采樣中心選取的影響[24]。為剔除采樣線上的建筑物點,對點云中任一點,計算統計高程大于該點的鄰域點數與全部鄰域點數的比值,比值越小,該點為地面點的可能性越大。對每條采樣線進行粗過濾后,采樣線上最左端點為采樣中心,距采樣中心最近的點為基準點。

1.3.2 基于斷裂線的采樣線點云遍歷

對任一條采樣線精細濾波,基于斷裂線約束,結合采樣中心點、基準點和各采樣點構建有向向量計算夾角[25],濾除地物點。確定采樣中心點和基準點后,從原始采樣線最左端開始遍歷,由近及遠按點云順序計算坡度ρ,得到采樣線上的地面點。坡度ρ的計算公式為

(11)

圖6為基于斷裂線的采樣線點云遍歷示意圖,具體的遍歷步驟如下。

圖6 采樣線點云遍歷Fig.6 A sampling line point cloud traversal

∥輸入:采樣中心點O、基準點B、采樣點數n

fori←0 ton

ifPi==斷裂點

保留Pi

ifPi!=斷裂點

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

ifρ≤ρthreshold(閾值)

保留Pi,Pi=B,i++

else

濾除Pi,i++

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

基于斷裂線約束的采樣線濾波結果如圖7所示,由圖7(b)可知,本文方法能準確區分采樣線上藍色地面點與紅色非地面點,確保后續處理順利執行。

圖7 采樣線濾波前后對比結果Fig.7 Comparison before and after sampling lines filtering

1.4 改進的二次曲面擬合

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

圖8 格網劃分Fig.8 Grids division

(12)

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

(13)

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

2 試驗結果與分析

為驗證本文方法的有效性,采用國際攝影測量與遙感學會(ISPRS)工作組提供的數據集和國內成都山區點云兩套數據進行測試。參考數據集是通過手動過濾點云數據集生成的,地形樣本中的每個點被分類為地面點或地物點,后續將手動分類結果作為真值檢驗方法精度。

2.1 ISPRS試驗數據

對于ISPRS提供的數據,本文選擇了4組具有不同地形特征的數據集對方法性能進行測試,詳細信息見表1。

表1 地形特征

2.1.1 試驗結果

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

表2 場景誤差分析

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

圖9 濾波結果Fig.9 Filtering results

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

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

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

2.1.2 參數設置

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

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

圖13 坡度閾值對應的總誤差Fig.13 Total errors for each slope threshold

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

圖14 格網邊長對應的總誤差Fig.14 Total errors for each grid size

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

圖15 曲率計算時鄰域范圍對應的總誤差Fig.15 Total error corresponding to different neighborhood in curvature calculation

2.1.3 消融研究

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

表3 添加斷裂線約束前后的誤差比較

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

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

表5 二次曲面擬合改進前后的誤差比較

2.2 成都山區試驗數據

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

圖16 成都山區機載點云數據Fig.16 Point clouds from the city of Chengdu

圖17 成都山區點云濾波結果Fig.17 The filtering results of the Chengdu data

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

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

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

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

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

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

為進一步探討斷裂線約束的重要性,本文分別計算圖18中A1、A2、A3區域參考DEM與有無斷裂線約束DEM之間的高差,定量對比有無斷裂線約束下的濾波后DEM精度,結果如圖20所示。誤差主要分布在斷裂線周圍,且對于Ⅰ類斷裂線和Ⅱ類斷裂線所在區域,無斷裂線約束的DEM與參考DEM的高差更大,表明基于斷裂線約束生成的DEM質量更高,在地形陡峭處表現更好。

3 結論與展望

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