代震,何榮,白偉森

(河南理工大學(xué) 測繪與國土信息工程學(xué)院,河南 焦作 454000)

0 引言

隨著實(shí)景三維中國建設(shè)的需求,對數(shù)字高程模型(digital elevation model,DEM)的精度要求越來越高,與以往為災(zāi)害監(jiān)測、地形繪圖等方向的支撐不同,精細(xì)DEM更能參與到智慧城市的建設(shè)中,關(guān)系到城市級實(shí)景三維模型的質(zhì)量。其中三維激光的廣泛應(yīng)用,使得獲取地面目標(biāo)的精確點(diǎn)云變得更加容易,點(diǎn)云處理技術(shù)已成為測繪領(lǐng)域的一個(gè)關(guān)鍵研究方向[1-3]。

現(xiàn)有的點(diǎn)云濾波方法一般分為3類。一是基于坡度的方法[4],該方法原理是以任意點(diǎn)構(gòu)建圓錐體,根據(jù)點(diǎn)云數(shù)據(jù)中地面坡度的大小分離地面點(diǎn)和非地面點(diǎn)。原理相對簡單,坡度濾波核運(yùn)算效率較高,但閾值設(shè)置困難,大多是依靠經(jīng)驗(yàn)量,同時(shí)算法在地形復(fù)雜的區(qū)域表現(xiàn)不佳。二是基于表面的方法,代表算法有漸進(jìn)加密三角網(wǎng)濾波算法(progressive tin densification,PTD)、布料模擬濾波算法(cloth simulation filter,CSF)等[5-9]。PTD算法原理是隨機(jī)選取局部區(qū)域最低點(diǎn)作為種子點(diǎn),然后根據(jù)上述種子點(diǎn)構(gòu)建稀疏三角網(wǎng)模型,通過判斷點(diǎn)云數(shù)據(jù)是否滿足距離和角度閾值來分離地面點(diǎn)和非地面點(diǎn)。該算法能夠較好地保留地形陡峭區(qū)域的地面點(diǎn),是目前較穩(wěn)健的算法之一。陳琳等[10]在2014年運(yùn)用高程統(tǒng)計(jì)的方法,進(jìn)一步升級基于先驗(yàn)經(jīng)驗(yàn)的PTD濾波算法,提高了濾波效率。凌曉春[11]在2020年引入薄板樣條曲線插值法,通過TPS中的彎曲能力增長值對PTD算法進(jìn)行改進(jìn),減少了濾波誤差。CSF算法原理是想象一個(gè)可以改變剛性系數(shù)的布料,緩緩下沉到倒置的地形上,以此得到地面點(diǎn)。該算法參數(shù)較少且易于設(shè)置,在平面地區(qū)具有極佳的濾波效果,但無法應(yīng)對高程變化劇烈的復(fù)雜地形。石壯等[12]在2022年通過虛擬格網(wǎng)劃分不同地形,分別采用對應(yīng)的參數(shù)進(jìn)行濾波,比單一CSF濾波更適用于混合地形。王佳雯等[13]提出一種基于地面點(diǎn)歸一化的濾波改進(jìn)方法,顯著提高濾波精度。三是基于形態(tài)學(xué)的方法。該方法通常采用數(shù)學(xué)應(yīng)用中的膨脹和腐蝕運(yùn)算,以固定的窗口分離地面點(diǎn)和非地面點(diǎn)。Zhang等[14]在2003年提出了一種經(jīng)典的漸進(jìn)式形態(tài)學(xué)濾波算法,該算法關(guān)系窗口尺寸、地形坡度、高差閾值等參數(shù),由于閾值的設(shè)置較為困難,使其一方面難以去除近地面點(diǎn),另一方面容易導(dǎo)致陡峭區(qū)域過度濾波。隋立春等[15]在數(shù)學(xué)形態(tài)學(xué)“開”算子的基礎(chǔ)上,首次提出“帶寬”的概念,通過增加該參數(shù)以提高最終濾波效果。苗啟廣等[16]在經(jīng)典形態(tài)學(xué)算法的基礎(chǔ)上,針對閾值設(shè)置問題,通過建立分塊區(qū)域預(yù)測地形坡度,可以根據(jù)區(qū)域地形起伏情況自適應(yīng)地調(diào)整閾值,得到最終結(jié)果。Pingel等[17]在前人研究的基礎(chǔ)上提出簡單形態(tài)學(xué)濾波算法(simple morphological filter,SMRF),利用輸入的最大濾波窗口尺寸、地面高程、地形坡度3個(gè)參數(shù)生成臨時(shí)地表,分類原始激光點(diǎn)云,算法更適用于城市地形,在非城市區(qū)域?yàn)V波誤差較大。

因此,針對上述SMRF算法在非城市區(qū)域?yàn)V波效果的不佳表現(xiàn),以及地面平面擬合算法在平緩地形的濾波優(yōu)勢,本文提出了一種基于地面平面擬合和簡單形態(tài)學(xué)濾波結(jié)合的點(diǎn)云濾波算法。算法能夠在抑制Ⅰ類誤差增大的前提下進(jìn)一步剔除近地面點(diǎn),大幅降低Ⅱ類誤差,提高生成的DEM質(zhì)量。

1 算法原理

從宏觀層面上來看,多數(shù)算法無法應(yīng)對大尺寸、多變化的地形地貌,尤其是在面對高程變化劇烈的陡坡區(qū)域,濾波誤差遠(yuǎn)遠(yuǎn)超出預(yù)料;反而對于較為精細(xì)的地物劃分,具有較好的濾波效果?；诖祟愃枷?在考慮到點(diǎn)云密度、曲率等多方面因素的影響權(quán)重,將劇烈變化的高程壓縮到較小范圍的輕微起伏也變成最有效策略。

本文具體算法流程如圖1所示,提出的濾波算法包括4個(gè)改進(jìn)步驟:簡單形態(tài)學(xué)粗濾波、DEM輔助的高程歸一化、地面平面擬合算法和空間向量后處理。

1.1 簡單形態(tài)學(xué)粗濾波

漸進(jìn)式形態(tài)學(xué)濾波算法是直接對原始點(diǎn)云進(jìn)行處理,分離出的非地面點(diǎn)主要依靠形態(tài)學(xué)開運(yùn)算的獨(dú)立篩選,處理過程不僅耗時(shí)耗力且容易造成誤判。而簡單形態(tài)學(xué)算法使用多尺寸窗口構(gòu)建柵格,以此生成臨時(shí)地形表面,通過多次迭代使表面吻合實(shí)際地形,最后統(tǒng)一進(jìn)行地面點(diǎn)與非地面點(diǎn)的劃分。

算法首先創(chuàng)建一個(gè)名為lastsurface的ZI副本。保留每個(gè)像元內(nèi)所有激光雷達(dá)點(diǎn)的最低高程,取代像元內(nèi)剩余點(diǎn)云的高程,判斷柵格單元內(nèi)點(diǎn)云數(shù)量為0時(shí)使用圖像修復(fù)技術(shù)插值。構(gòu)建圓盤形結(jié)構(gòu),對lastsurface進(jìn)行數(shù)學(xué)形態(tài)學(xué)開運(yùn)算,如果原數(shù)據(jù)高程與開運(yùn)算結(jié)果高程之差大于設(shè)定參數(shù)(機(jī)載激光數(shù)據(jù)默認(rèn)為0.3),把該點(diǎn)當(dāng)做空單元,使用圖像修復(fù)技術(shù)進(jìn)行插值,小于設(shè)定值,就把該點(diǎn)看做地面單元。激光數(shù)據(jù)需要不斷迭代,迭代步長為一個(gè)單元,需要提前設(shè)置好最大尺寸。處理后創(chuàng)建臨時(shí)DEM,保留ZI中被確定為地面的單元,然后根據(jù)原始點(diǎn)與臨時(shí)DEM間的高差進(jìn)行分類。

1.2 DEM輔助的高程歸一化

在實(shí)際實(shí)驗(yàn)中,根據(jù)最大最小點(diǎn)云坐標(biāo)的高程歸一化結(jié)果誤差較大,容易受到植被和建筑物的影響。因此采用SMRF濾波生成的粗DEM進(jìn)行高程歸一化,歸一化原理如圖2所示。

首先采用不規(guī)則三角網(wǎng)生成粗DEM,然后統(tǒng)計(jì)原始點(diǎn)云坐標(biāo),構(gòu)建格網(wǎng)映像,通過待處理點(diǎn)云坐標(biāo)投影格網(wǎng),找到對應(yīng)的DEM擬合高程數(shù)據(jù),與待處理點(diǎn)云的原始高程作差,得到新的歸一化值。同時(shí),保留激光腳點(diǎn)的原始高程值,便于處理后反歸一化。

1.3 地面平面擬合算法

高程歸一化后的點(diǎn)云趨近于平緩地形,這些區(qū)域的地形可以看作連續(xù)分布的曲面,在相對較小的區(qū)域內(nèi),該曲面可以近似為等效平面。本文在地面平面擬合的基礎(chǔ)上,加入漸進(jìn)移動(dòng)窗口和格網(wǎng)的改進(jìn)步驟,同時(shí)以每個(gè)格網(wǎng)內(nèi)的最低點(diǎn)作為數(shù)據(jù)集擬合平面模型,以此形成多個(gè)平面切片。算法步驟如下。

步驟1:建立索引機(jī)制,以快速查詢每個(gè)點(diǎn)所在的網(wǎng)格和每個(gè)網(wǎng)格包含的點(diǎn)。點(diǎn)與網(wǎng)格間的索引關(guān)系的計(jì)算如式(1)所示。

(1)

式中:(X,Y)為網(wǎng)格號(hào);(x,y)為點(diǎn)云的平面坐標(biāo);(xmin,ymin)為整個(gè)數(shù)據(jù)集的最小平面坐標(biāo);m為網(wǎng)格單元,即移動(dòng)窗口大小;INT表示對計(jì)算結(jié)果向下取整。

步驟2:種子點(diǎn)選取。在每個(gè)格網(wǎng)內(nèi),索引最低點(diǎn)構(gòu)建種子數(shù)據(jù)集S:si{i=1,2,…,n}(n≥3),定義平面為

ax+by+cz+d=0

(2)

假設(shè)w=(abc)T,x=(xyz)T,式(2)簡化為式(3)。

-wTx=d

(3)

步驟3:擬合平面模型。通過初始種子點(diǎn)集S構(gòu)建協(xié)方差矩陣C∈R3×3,如式(4)、式(5)所示。

(4)

Cx=λx

(5)

(6)

單一閾值的設(shè)定具有局限性,窗口變換的同時(shí)漸進(jìn)縮小地面點(diǎn)提取閾值D0,濾除與當(dāng)前模型距離較大的點(diǎn)。通過大量實(shí)驗(yàn),設(shè)置D0初始值為5、步長為0.8時(shí),提取效果最佳。平面距離D與閾值D0比較,若D≤D0,將P點(diǎn)作為地面點(diǎn)合并到此面集;若D≥D0,則作為非地面點(diǎn)濾除,重復(fù)數(shù)次后得到較為精確的平面點(diǎn)集。

1.4 空間向量后處理

在實(shí)際應(yīng)用中,平面模型的構(gòu)建容易受到其他地物的影響。為了提高地平面擬合的精度,利用區(qū)域空間向量投影,驗(yàn)證所擬合地平面的準(zhǔn)確性,避免將一些異常平面(如建筑頂面、立面等)分割成地面。首先,統(tǒng)計(jì)各面片集合點(diǎn)云數(shù)量,從多到少依次排列,以點(diǎn)云數(shù)量最多的面集作為標(biāo)準(zhǔn)平面,設(shè)置標(biāo)準(zhǔn)法向量q;然后,對后續(xù)擬合出的地平面集合進(jìn)行判斷,空間向量二面角公式為式(7)。

(7)

式中:p為待測平面集合的法向量。

由此計(jì)算出二面角,通過與閾值比較,若小于閾值,則與標(biāo)準(zhǔn)平面合并。針對二面角閾值的設(shè)置問題,太大容易造成多個(gè)平面過度合并;設(shè)置太小容易造成地平面過度分割。提出一種根據(jù)各面片集合點(diǎn)云密度判斷的解決方法,待測面片集合的點(diǎn)云密度越接近標(biāo)準(zhǔn)平面點(diǎn)云密度,越可能是真實(shí)地平面,公式為式(8)。

(8)

式中:ρ為標(biāo)準(zhǔn)平面點(diǎn)云密度;ρi為待測面片點(diǎn)云密度;α0為二面角閾值初值;αi為待測平面二面角確定閾值。

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

2.1 數(shù)據(jù)來源

實(shí)驗(yàn)數(shù)據(jù)分為兩類。第一類實(shí)驗(yàn)數(shù)據(jù)采用飛馬D-LiDAR2000系統(tǒng)于2021年4月14日獲取,采集地點(diǎn)位于河南理工大學(xué)南校區(qū),激光數(shù)據(jù)如圖3所示,點(diǎn)云密度為7.49/m2,同時(shí)采集實(shí)驗(yàn)區(qū)域GPS點(diǎn)位進(jìn)行算法驗(yàn)證。第二類實(shí)驗(yàn)所用數(shù)據(jù)來自國際攝影測量與遙感協(xié)會(huì)(ISPRS)官方網(wǎng)站發(fā)布的標(biāo)準(zhǔn)點(diǎn)云數(shù)據(jù),數(shù)據(jù)采集地點(diǎn)位于 Vaihingen/Enz測試場和Stuttgart市中心,選取8組具有代表性的測試數(shù)據(jù)進(jìn)行實(shí)驗(yàn),如表1所示。兩組實(shí)驗(yàn)數(shù)據(jù)包含建筑物、公路、植被、陡坡等不同地物特征,基本對應(yīng)常見地形,數(shù)據(jù)集已手動(dòng)分離出地面點(diǎn)和非地面點(diǎn)。

表1 研究區(qū)域點(diǎn)云基本情況

圖3 河南理工大學(xué)南校區(qū)

2.2 理工校區(qū)DEM精度評價(jià)

本文采用不規(guī)則三角網(wǎng)方法處理濾波后的實(shí)驗(yàn)數(shù)據(jù),以SMRF算法和改進(jìn)算法提取的地面點(diǎn)分別生成DEM模型,由兩種算法對應(yīng)的DEM精度反映二者的濾波精度。已知41個(gè)GPS地面檢查點(diǎn)的三維坐標(biāo),將其視為真實(shí)地面高程,從生成的DEM模型中提取相同點(diǎn)位的擬合高程值。取擬合高程值與真實(shí)高程值差值的中誤差(root mean square error,RMSE)和平均絕對誤差(mean absolute error,MAE)作為濾波精度的判別依據(jù),并分別與實(shí)際地形進(jìn)行擬合。

由圖4中殘差分布可知,兩種濾波算法濾波后DEM的殘差分布趨勢大體相似,但是SMRF算法對應(yīng)的DEM的殘差比改進(jìn)算法波動(dòng)振幅大,改進(jìn)算法的波動(dòng)更為平緩,產(chǎn)生的誤差更小。

圖4 殘差分布

由表2可知,與SMRF算法相比,改進(jìn)算法的中誤差降低了3.3 cm,平均絕對誤差降低了2.0 cm,證明改進(jìn)算法能夠有效地降低點(diǎn)云誤差。將兩種算法分別與實(shí)際地形線性擬合,如圖5、圖6所示,改進(jìn)算法生成的DEM模型與實(shí)際地形的擬合效果更好、相關(guān)度更強(qiáng),更加逼近于真實(shí)地形。

表2 DEM指標(biāo)統(tǒng)計(jì) m

圖5 SMRF算法與實(shí)際地形線性擬合

圖6 改進(jìn)算法與實(shí)際地形線性擬合

2.3 國際標(biāo)準(zhǔn)數(shù)據(jù)評價(jià)

以國際標(biāo)準(zhǔn)數(shù)據(jù)進(jìn)行多種算法比較,本文采用2003年美國國際攝影與遙感協(xié)會(huì)提出由混淆矩陣推導(dǎo)的一種交叉表評價(jià)體系,如表3所示。其中,Ⅰ類誤差(TⅠ)是指地面點(diǎn)中被錯(cuò)誤地劃分到非地考慮到CSF濾波算法誤差較為典型,分別使用CSF算法、SMRF算法和本文算法對8組數(shù)據(jù)進(jìn)行實(shí)驗(yàn)并比較,得到的Ⅰ類誤差、Ⅱ類誤差和總誤差結(jié)果如表4所示。根據(jù)表4可以看出,本文算法的Ⅰ類誤差、Ⅱ類誤差和總誤差分別為2.54%、7.47%和3.06%,均小于另兩種算法。為了更好地驗(yàn)證實(shí)驗(yàn)效果,這里展示了地形顯著、效果明顯的Samp11、Samp23、Samp52濾波結(jié)果,并進(jìn)行對比分析。

表4 各研究區(qū)濾波結(jié)果統(tǒng)計(jì) %

面點(diǎn)的數(shù)量占地面點(diǎn)總數(shù)的百分比,Ⅱ類誤差(TⅡ)表示非地面點(diǎn)中被錯(cuò)誤地劃分到地面點(diǎn)的數(shù)量占非地面點(diǎn)總數(shù)的百分比,總誤差(TⅢ)是分類結(jié)果與參考數(shù)據(jù)不一致的概率。

進(jìn)一步分析,本文算法在各個(gè)實(shí)驗(yàn)區(qū)的Ⅰ類誤差均小于8.1%,表明該算法在處理陡坡、大型建筑物、植被覆蓋和數(shù)據(jù)間斷地形環(huán)境時(shí)容易得到地面點(diǎn),同時(shí)能夠保留更多的地形細(xì)節(jié);處理實(shí)驗(yàn)區(qū)Samp51、Samp52時(shí)Ⅱ類誤差較大,這說明本文算法在處理低矮植被和河流河岸區(qū)域時(shí)容易將非地面點(diǎn)判別為地面點(diǎn),而Ⅱ類誤差通常較容易通過人工剔除。

Samp11數(shù)據(jù)是建立在斜坡上的城市區(qū)域。如圖7(a)的圓圈標(biāo)注所示,在斜坡坡度變化最大區(qū)域和較高建筑區(qū)域,CSF算法將地面點(diǎn)云當(dāng)成非地面點(diǎn)剔除,出現(xiàn)點(diǎn)云空洞較多,原因是密集城區(qū)的集聚,布料硬度參數(shù)和迭代次數(shù)都設(shè)置的較小,對于斜坡區(qū)域來說,地形起伏和建筑落差較大,地面點(diǎn)容易被當(dāng)作非地面點(diǎn)剔除。如圖7(b)的圓圈標(biāo)注所示,SMRF算法將一些低矮建筑物和植被區(qū)域附近的地面點(diǎn)錯(cuò)誤剔除,降低濾波精度,原因在于低矮地物過于貼合地面,干擾濾波過程。而本文方法在相應(yīng)坡度較大的區(qū)域,仍保留地形細(xì)節(jié)。

圖7 不同算法對Samp11的濾波效果對比

Samp23數(shù)據(jù)是典型的大型復(fù)雜建筑區(qū)。如圖8(a)所示,右上角圓圈標(biāo)注是明顯的低矮建筑物,CSF算法錯(cuò)誤地把非地面點(diǎn)判斷成地面點(diǎn);而左下角圓圈標(biāo)注,由于建筑臺(tái)階依托的地勢較低,遠(yuǎn)低于正常路面,形成坡度變化較大的地形條件,CSF算法濾波有明顯的點(diǎn)云空洞出現(xiàn)。如圖8(b)的圓圈標(biāo)注所示,SMRF算法錯(cuò)誤地將地面點(diǎn)判斷成非地面點(diǎn),該區(qū)域地勢略高于路面。對比結(jié)果顯示本文算法在大型建筑區(qū)具有很好的濾波效果,對略低于或高于正常地面的建筑臺(tái)階也能很好保留,并作為地形地物特征。

圖8 不同算法對Samp23的濾波效果對比

Samp52數(shù)據(jù)是流經(jīng)河流的山坡地帶。如圖9所示,圓圈標(biāo)注區(qū)域都在坡峰附近,其側(cè)面過于陡峭,地形變化劇烈,CSF算法和SMRF算法在應(yīng)對該區(qū)域時(shí)都容易把地面點(diǎn)錯(cuò)誤判斷成非地面點(diǎn);同時(shí)河岸附近存在大量低矮植被和建筑物,也對濾波過程造成了干擾。對比結(jié)果顯示,雖然本文算法在該區(qū)域具有相對更好的濾波效果,但也受到陡峭地形和低矮地物的干擾,出現(xiàn)較大的Ⅱ類誤差。

圖9 不同算法對Samp52的濾波效果對比

為了更客觀、清楚地評判改進(jìn)算法的優(yōu)劣性,對本文算法的濾波結(jié)果與已測試過的其他4種濾波算法的精度進(jìn)行比較分析,包括CSF濾波算法、移動(dòng)曲面濾波算法、PTD濾波算法和SMRF濾波算法。如圖10所示,在5種濾波算法對8個(gè)不同地形的實(shí)驗(yàn)中,CSF算法在Samp12、Samp23、Samp31、Samp51實(shí)驗(yàn)區(qū)總誤差小于11.5%,Samp53區(qū)域總誤差為6.6%;移動(dòng)曲面算法在Samp12、Samp31、Samp42研究區(qū)總誤差小于10.1%;PTD算法在Samp31、Samp42區(qū)域總誤差小于9%,在實(shí)驗(yàn)區(qū)Samp51、Samp52、Samp53總誤差小于5.2%;SMRF算法在區(qū)域Samp12、Samp23、Samp31、Samp52的總誤差小于11.2%;本文算法在實(shí)驗(yàn)區(qū)總誤差均小于5.4%,在實(shí)驗(yàn)區(qū)Samp12、實(shí)驗(yàn)區(qū)Samp42的總誤差接近于0。比較發(fā)現(xiàn),本文方法的平均總誤差最小,同時(shí)各個(gè)實(shí)驗(yàn)區(qū)的Ⅰ類誤差始終小于8.1%,證明該算法具有很好的適用性,能夠有效地降低點(diǎn)云誤差。

圖10 5種算法的濾波精度比較

3 結(jié)束語

本文結(jié)合簡單形態(tài)學(xué)濾波,改進(jìn)了地面平面擬合算法;利用地面平面擬合算法對歸一化點(diǎn)云的濾波優(yōu)勢,顯著提高了算法精度;通過漸進(jìn)移動(dòng)窗口構(gòu)建多個(gè)平面模型,減少迭代步驟,克服單一閾值的局限性;針對多個(gè)平面的合并問題,采用空間向量的后處理方法,以點(diǎn)云密度控制二面角閾值,避免過度分割。

算法濾波精度得到有效提高。與SMRF算法比較,本文算法的中誤差降低了3.3 cm,平均絕對誤差降低了2.0 cm。與CSF算法、移動(dòng)曲面算法、PTD算法和SMRF算法比較,本文算法的平均總誤差分別減少8.8%、12.5%、6.5%和9.2%,能夠在抑制Ⅰ類誤差增大的同時(shí)大幅降低Ⅱ類誤差,證明該算法能夠有效地降低點(diǎn)云誤差。

濾波方法適用于多數(shù)地形條件。實(shí)驗(yàn)區(qū)包含建筑物、公路、植被、陡坡等不同地物特征,基本對應(yīng)常見地形,本文算法在實(shí)驗(yàn)區(qū)總誤差均小于5.4%,在實(shí)驗(yàn)區(qū)Samp12、實(shí)驗(yàn)區(qū)Samp42的總誤差接近于0,同時(shí)各個(gè)實(shí)驗(yàn)區(qū)Ⅰ類誤差始終小于8.1%,表明本文算法能夠適用于多數(shù)地形條件。