羅為東,甘 淑,袁希平,高 莎,畢 瑞,袁新悅

(1.昆明理工大學(xué)國土資源工程學(xué)院,云南 昆明 650093；2.云南省高校高原山地空間信息測繪技術(shù)應(yīng)用工程研究中心,云南 昆明 650093；3.滇西應(yīng)用技術(shù)大學(xué)地球科學(xué)與工程技術(shù)學(xué)院,云南 大理671006)

1 引 言

近年來隨著科學(xué)技術(shù)的發(fā)展,復(fù)雜山區(qū)地形地貌提取分析朝著精細(xì)化方向發(fā)展,如何高效、快速獲取高精度地形地貌空間數(shù)據(jù)成為了分析地形地貌特征的關(guān)鍵因素。無人機(jī)航測系統(tǒng)在對生產(chǎn)地面表面模型(Digital Surface Model,DSM)、數(shù)字正射影像圖(Digital Orthophoto Map,DOM)、構(gòu)建數(shù)字高程模型(Digital Elevation Model,DEM)等具有快速、高效的特點(diǎn)。地形數(shù)學(xué)產(chǎn)品均是對地形信息的詳細(xì)表達(dá)。其中數(shù)字高程模型的快速獲取,并滿足相應(yīng)的精度要求,特別是在困難山區(qū)對數(shù)字高程模型的處理,是山區(qū)數(shù)字化發(fā)展的重要趨勢。目前,國內(nèi)外學(xué)者利用無人機(jī)遙感技術(shù)采集構(gòu)建DEM對地形地貌等空間特征分析展開了研究。不難發(fā)現(xiàn),構(gòu)建DEM對于無人機(jī)航測技術(shù),主要是要對密集影像匹配(Dense image matching,DIM)點(diǎn)云進(jìn)行濾波處理,對此GordanaJakovljevic[1]等通過LiDAR點(diǎn)云和密集匹配點(diǎn)云進(jìn)行濾波處理,構(gòu)建DEM用于降低洪水危險,并計(jì)算DEM的RMSE,結(jié)果表明平原地區(qū)適宜應(yīng)用LiDAR采集點(diǎn)云,而丘陵和山區(qū)密集匹配點(diǎn)云構(gòu)建DEM的RMSE低于LiDAR點(diǎn)云。Mustafa Zeybek[2]等應(yīng)用四種濾波算法結(jié)果提取地形特征,結(jié)果表明粗糙度和復(fù)雜地貌提取四種濾波算法表現(xiàn)出相似性,CSF濾波算法在平坦地形區(qū)分地面點(diǎn)和非地面點(diǎn)正確率達(dá)到93 %。葉立志[3]針對DIM點(diǎn)云本身的復(fù)雜性,提出了一種基于交叉線元分割的密集匹配點(diǎn)云濾波算法,結(jié)果表明濾波結(jié)果構(gòu)建DEM中誤差接近1m,相對于商業(yè)軟件Terrasolid的濾波結(jié)果和濾波效率更優(yōu)。JiXian[4]等應(yīng)用ISPRS第三個工作組的七個數(shù)據(jù)集,首先構(gòu)建完初試TIN后,再SUSC擴(kuò)展得到更多地面點(diǎn),最后進(jìn)行對TIN進(jìn)行迭代加密,結(jié)果表明與PTD相比,所提出的方法能夠預(yù)景觀的不連續(xù)性,并將遺漏誤差和總誤差分別減少了大約10 %和6 %,降低了后期校正結(jié)果所需的人工操作成本。ShengNie[5]等對PTD濾波算法中建立TIN進(jìn)行改進(jìn)并改變原有的迭代判斷標(biāo)準(zhǔn),在ISPRS提供LiDAR點(diǎn)云數(shù)據(jù)集進(jìn)行處理,結(jié)果表明改進(jìn)的PTD方法能夠把I類誤差、Ⅱ誤差和總誤差分別減少10.26 %、0.79 %和8.07 %。為過濾機(jī)載LiDAR點(diǎn)云離散數(shù)據(jù)集濾波供了決方案。高廣[6]等針對山區(qū)LiDAR點(diǎn)云的特點(diǎn),對PTD濾波首先應(yīng)用隨機(jī)格網(wǎng)搜索算法獲取更多精準(zhǔn)初試種子點(diǎn),其次利用地形預(yù)測角的判斷準(zhǔn)則提高地線附近濾波精度,最后在有大量地形斷裂線山區(qū)進(jìn)行驗(yàn)證,結(jié)果表明改進(jìn)PTD濾波算法有效地保留了山區(qū)的地形斷裂線特征。

對于祿豐恐龍谷南緣山區(qū)特殊地形地貌,采用低空無人機(jī)測量技術(shù)獲取的影像數(shù)據(jù),并構(gòu)建DIM點(diǎn)云,首先考慮山頂、山體兩側(cè)和地面存在高程差,選擇漸進(jìn)加密三角網(wǎng)濾波算法進(jìn)行處理。其次針對實(shí)驗(yàn)區(qū)特殊地理位置和地貌形態(tài),改進(jìn)傳統(tǒng)漸進(jìn)加密三角網(wǎng)濾波算法,最后對兩種濾波算法處理結(jié)果進(jìn)行直觀分析。較為完整得到實(shí)驗(yàn)區(qū)地面點(diǎn)數(shù)據(jù)可為后續(xù)解譯地貌特征信息和提取地形特征信息奠定基礎(chǔ)。

2 PTD改進(jìn)方法研究

2.1 DIM點(diǎn)云PTD濾波方法

目前無人機(jī)航測制作的主要數(shù)據(jù)之一是點(diǎn)云數(shù)據(jù),以點(diǎn)集的形式表達(dá)被測物體的屬性、位置信息和紋理信息。影像密集匹配獲取點(diǎn)云,即DIM點(diǎn)云,不可透過地表的地物獲取地面點(diǎn),獲取的地表點(diǎn)云與建筑物、信號塔和植被等地物構(gòu)成地面表面模型(Digital Surface Model,DSM),存在噪聲點(diǎn)及高程異常點(diǎn)。密集匹配點(diǎn)云濾波,在點(diǎn)云分類、三維建模、單體化提取等方面有重要作用。國內(nèi)外研究者[7-14]在點(diǎn)云濾波方面開展了大量工作,依據(jù)主要濾波原理形成了基于坡度、基于形態(tài)學(xué)、基于曲面擬合、基于分割、基于機(jī)器學(xué)習(xí)、布料模擬濾波和漸進(jìn)加密不規(guī)則三角網(wǎng)算法等主要濾波方法?；诓灰?guī)則三角網(wǎng)的濾波算法,在丘陵地區(qū)、林區(qū)和山區(qū)等對點(diǎn)云數(shù)據(jù)有較好的適用性。經(jīng)典的漸進(jìn)加密不規(guī)則三角網(wǎng)(Progressive TIN Densification,PTD)濾波算法獲取實(shí)驗(yàn)區(qū)地面點(diǎn)。PTD算法在商用軟件Terrasolid的TerraScan模塊中采用,對實(shí)驗(yàn)區(qū)濾波主要進(jìn)行四步:①把實(shí)驗(yàn)區(qū)數(shù)據(jù)劃分規(guī)則格網(wǎng),認(rèn)定格網(wǎng)最低點(diǎn)為地面點(diǎn),在實(shí)驗(yàn)區(qū)選取初試種子點(diǎn)；②利用種子點(diǎn)構(gòu)建不規(guī)則三角網(wǎng)(Triangulated Irregular Network,TIN)；③判斷地面點(diǎn),如圖1所示,P點(diǎn)為非種子點(diǎn),底面的三角平面頂點(diǎn)分別是V1、V2、V3,d表示非種子點(diǎn)到三角面的距離,α1、α2、α3是非種子點(diǎn)與底面三角形頂點(diǎn)連線的夾角,通過計(jì)算點(diǎn)P到三角形的距離d及α1、α2、α3的最大值分別與距離和角度的閾值比較,若小于閾值則定為地面點(diǎn),歸入三角網(wǎng)中。④迭代循環(huán)步驟②、③直至分類完所有實(shí)驗(yàn)區(qū)點(diǎn)云。

圖1 判斷地面點(diǎn)

2.2 PTD算法的改進(jìn)技術(shù)

PTD濾波算法每次迭代過程中,都對其余各點(diǎn)到所在三角形的反復(fù)角和反復(fù)距離進(jìn)行閾值判斷,將滿足條件的點(diǎn)加到TIN,初始構(gòu)建的TIN對后續(xù)濾波判斷影響較大。PTD濾波算法把尺寸大于研究區(qū)最大建筑物尺寸的格網(wǎng)中的最低點(diǎn)視為地面點(diǎn),從該角度出發(fā)又考慮實(shí)驗(yàn)區(qū)山體兩側(cè)溝壑叢生,且對于山地濾波,山體與地面存在高程差,山體的地面點(diǎn)更易被誤識別為非地面點(diǎn)。因此在PTD濾波算法選取種子點(diǎn)時,加入山脊線和山谷線交匯的地形特征點(diǎn),通過地形特征點(diǎn)使山體兩側(cè)構(gòu)建密集且符合實(shí)際地形的三角網(wǎng),使山體兩側(cè)的點(diǎn)易通過角度和距離的閾值。針對實(shí)驗(yàn)區(qū)特殊地形分布情況,改進(jìn)種子點(diǎn)選擇方法,使點(diǎn)位分類在構(gòu)建TIN中符合實(shí)際情況。如2圖所示,C1～C5是初始種子點(diǎn),A和B是在脊-谷的位置選擇特征點(diǎn),AB連線是構(gòu)建三角網(wǎng)底邊,點(diǎn)G1和G2是三角網(wǎng)下部的點(diǎn),與B點(diǎn)在同一三角網(wǎng)但高程值低于A點(diǎn)。G2點(diǎn)為非地面點(diǎn),G1為地面點(diǎn),所以仍不能將三角網(wǎng)的下部點(diǎn)均視為地面點(diǎn)。

為此需要不斷剔除最遠(yuǎn)點(diǎn)和增加點(diǎn)構(gòu)建移動曲面,確保點(diǎn)云規(guī)則傳遞并覆蓋整個實(shí)驗(yàn)區(qū),如3圖示,主要進(jìn)行:①把實(shí)驗(yàn)區(qū)固定格網(wǎng)按建筑物最大尺寸劃分；②確定實(shí)驗(yàn)區(qū)最外圍格網(wǎng)高程最低點(diǎn),以此為地面種子點(diǎn)；③點(diǎn)A、B、C、D是研究區(qū)四個角格網(wǎng)中心,以此構(gòu)網(wǎng),域內(nèi)全部點(diǎn)為中心建立移動格網(wǎng)；④如圖4示,進(jìn)行多尺度迭代把移動格網(wǎng)內(nèi)最低高程點(diǎn)作為地面種子點(diǎn)；⑤剔除重復(fù)點(diǎn),并和格網(wǎng)外地面種子點(diǎn)建立實(shí)驗(yàn)區(qū)地面種子點(diǎn)。

圖2 表示待定點(diǎn)識別

圖3 移動格網(wǎng)構(gòu)建

圖4 格網(wǎng)多尺度迭代

為避免格網(wǎng)缺失產(chǎn)生的錯誤分類,格網(wǎng)ABCD外部不再構(gòu)建移動格網(wǎng)。每個移動格網(wǎng)均依照實(shí)驗(yàn)區(qū)建筑物最大尺寸劃分,既保證了所有格網(wǎng)均有地面點(diǎn),又確保獲取地面種子點(diǎn)無需進(jìn)一步優(yōu)化,直接進(jìn)入三角網(wǎng)迭代判斷的步驟。

3 試驗(yàn)區(qū)及試驗(yàn)數(shù)據(jù)采集

3.1 試驗(yàn)區(qū)

研究區(qū)域位于云南祿豐恐龍國家地質(zhì)公園南緣山區(qū),隸屬云南楚雄彝族自治州,地理坐標(biāo)為N 24°51′33″～25°30′45″,E 101°38′06″～102°24′34″。測區(qū)范圍內(nèi)地貌類型復(fù)雜多樣,以構(gòu)造侵蝕地貌、方山地貌為主,地勢東北高,西南低,最高海拔2200 m,最低海拔1302 m,干濕分明。由于測區(qū)內(nèi)存在小型中生代紅色沉積盆地,成土母巖由中生代紫色砂頁巖和元古代的碳酸鹽巖交錯分布,受母巖影響,土壤類型呈帶分布,紫色和紅色土壤分布廣泛,測區(qū)位置如圖5所示。

圖5 實(shí)驗(yàn)區(qū)概況

3.2 試驗(yàn)數(shù)據(jù)采集

試驗(yàn)區(qū)影像數(shù)據(jù)獲取是通過DJI Phantom 4 RTK進(jìn)行影像數(shù)據(jù)采集。首先根據(jù)測區(qū)實(shí)際周邊地理環(huán)境,交通狀況,結(jié)合Google earth平臺對測區(qū)進(jìn)行航線規(guī)劃,其次通過前期航線規(guī)劃以及后期數(shù)據(jù)質(zhì)量要求,選擇合適的參數(shù)設(shè)置。鑒于本次數(shù)據(jù)獲取主要面向微地貌特征分析應(yīng)用,故本次飛行參數(shù)設(shè)置為:航向重疊度80 %,旁向重疊度80 %,平均飛行高度150 m。本次航測天氣條件良好,共采集392張影像,影像平均分辨率為0.07 m。

4 實(shí)驗(yàn)分析

4.1 基于DIM點(diǎn)云可視化分析

通過對實(shí)測影像進(jìn)行提取同名像點(diǎn)、特征匹配、計(jì)算旋轉(zhuǎn)矩陣、對選擇矩陣進(jìn)行畸變校正和光束法平差迭代完成影像密集匹配獲取DIM點(diǎn)云。點(diǎn)云密度為154.66個/m2,標(biāo)準(zhǔn)差為16.13 m。因無人機(jī)航測影像僅包含可見光波段,不具備透過地表已有地物如:植被、建筑物、信號塔等,即DIM點(diǎn)云含有噪聲點(diǎn)和高程異常點(diǎn),初步獲得地表DIM點(diǎn)云生成DSM,如圖6所示。

圖6 數(shù)字表面模型(DSM)

由圖6可知,恐龍谷南緣山區(qū)受亞熱帶低緯度高原季風(fēng)氣候影響山體溝壑交錯,常年干濕分明,也導(dǎo)致山體右側(cè)坡度通過目視解譯判斷接近90°。DSM中明顯高程異常有T1、T2和T3的三個區(qū)域,在數(shù)字正射影像(Digital Orthophoto Map:DOM)圖7中選定同樣三個區(qū)域,DOM包含真實(shí)色彩信息,可判斷出T1區(qū)域表示實(shí)測區(qū)域植被覆蓋區(qū)域,主要是低矮獨(dú)立樹和灌木。T2區(qū)域是山頂?shù)男盘査?信號塔在一定程度確保了免像控?zé)o人機(jī)對山區(qū)進(jìn)行航測時獲取POS的精度,但因信號塔自身高程值,后續(xù)構(gòu)建DEM同樣需先進(jìn)行點(diǎn)云濾波處理。T3區(qū)域是山腳的建筑物和蔬菜大棚,經(jīng)實(shí)地勘察該區(qū)域建筑物不超兩層,除上述三個區(qū)域的地表地物,實(shí)驗(yàn)區(qū)西北角和東南角分布的低矮灌木和農(nóng)作物均需要濾波處理。

圖7 正射影像圖(DOM)

4.2 基于PTD改進(jìn)的點(diǎn)云濾波處理成效對比分析

4.2.1 PTD濾波處理及去噪成效分析

實(shí)驗(yàn)通過“試錯法”調(diào)整閾值,選定最大地形坡度角為88°,迭代角度閾值設(shè)置為6°,迭代距離設(shè)置為1.4 m。70個初試種子點(diǎn)選取如圖8(a)示,種子點(diǎn)除了邊界點(diǎn),其余均勻分布于實(shí)驗(yàn)區(qū),由70個種子點(diǎn)構(gòu)建TIN如圖8(b)示,構(gòu)建三角網(wǎng)較為規(guī)整。

(a) 濾波種子點(diǎn)

(b) 種子點(diǎn)構(gòu)建不規(guī)則三角網(wǎng)

實(shí)驗(yàn)區(qū)采集DIM點(diǎn)云有17271784個點(diǎn),由初試種子點(diǎn)構(gòu)建TIN進(jìn)行PTD濾波,剔除非地面點(diǎn)并進(jìn)行冗余數(shù)據(jù)抽希,最后保留地面點(diǎn)5723274個點(diǎn)。

4.2.2 改進(jìn)PTD濾波處理技術(shù)應(yīng)用

對于上述PTD濾波算法對實(shí)驗(yàn)區(qū)造成山體部分區(qū)域地面點(diǎn)被誤識別為非地面點(diǎn),實(shí)驗(yàn)區(qū)位于滇中高原,山體頂部兩側(cè)與山底存在高程差,山頂和山脊區(qū)域易誤選為非地面點(diǎn)。因此,在PTD濾波算法基礎(chǔ)上選取種子點(diǎn)時,依據(jù)實(shí)驗(yàn)區(qū)特性,在ArcGIS中基于地表水流分析和幾何分析相結(jié)合的方法提取山脊線和山谷線,二者交匯的脊-谷地形特征點(diǎn)為改進(jìn)PTD濾波算法的種子點(diǎn)。圖9為主要山脊線、山谷線分布和脊-谷交匯點(diǎn)示意。

圖9 主要脊-谷線示意

對于地形特征點(diǎn)的選取精度評價需構(gòu)建一個評價模型。評價模型的構(gòu)建主要包含:①剔除冗余數(shù)據(jù),保留主要脊-谷特征線；②統(tǒng)計(jì)地形特征點(diǎn)數(shù)量,用M表示；③主要脊-谷特征線以1 m、5 m、10 m和15 m為線性單位構(gòu)建緩沖區(qū),計(jì)算落入緩沖區(qū)范圍內(nèi)地形特征點(diǎn)數(shù)量,壓蓋邊緣的點(diǎn)同樣算入,用Ni(i=1,5,10,15)表示其數(shù)量；④Qi(i=1,5,10,15)表示在緩沖區(qū)內(nèi)地形特征點(diǎn)數(shù)量占總點(diǎn)數(shù)的比,具體如式1,以此為地形特征點(diǎn)選取精度標(biāo)準(zhǔn)。

(1)

統(tǒng)計(jì)分析如表1所示,本次共提取486個地形特征點(diǎn),以1 m為距離的緩沖區(qū)中涵蓋27個點(diǎn),僅占總數(shù)的5 %,5 m為距離的緩沖區(qū)包含132個點(diǎn),占總數(shù)的27.16 %,15 m為距離的緩沖區(qū)有429個點(diǎn),占88.27 %。后續(xù)試驗(yàn)較少的種子點(diǎn)會造成濾波精度低,較多的種子點(diǎn)在處理過程中導(dǎo)致處理效率不高,故以10 m為緩沖距離覆蓋地形特征點(diǎn)為改進(jìn)PTD濾波算法的種子點(diǎn)。

表1 地形特征點(diǎn)提取精度評價

以圖8(a)初始種子點(diǎn)并疊加提取的脊-谷交匯的地形特征點(diǎn),二者共有300個點(diǎn)為種子點(diǎn),如圖10所示。通過此類改進(jìn)使選取種子點(diǎn)貼合實(shí)驗(yàn)區(qū)地形分布,并對選取種子點(diǎn)完成TIN構(gòu)建如圖10(b)所示。通過改進(jìn)方法構(gòu)建三角網(wǎng),位于山體兩側(cè)相對于初始構(gòu)建,三角網(wǎng)更為密集。改進(jìn)濾波算法并進(jìn)行冗余點(diǎn)抽希保留6692859個地面點(diǎn),相對于PTD濾波算法保留地面點(diǎn)更充足。

圖10 (a)地形特征點(diǎn)

圖10 (b)特征點(diǎn)構(gòu)建不規(guī)則三角網(wǎng)

4.2.3 改進(jìn)PTD去噪結(jié)果對比分析

對兩種濾波方法保留的地面點(diǎn)進(jìn)行目視解譯,對PTD濾波算法得到地面點(diǎn)結(jié)果和脊-谷交匯特征點(diǎn)疊加初始選取種子點(diǎn)結(jié)合構(gòu)建不規(guī)則三角網(wǎng)點(diǎn)云濾波結(jié)果均顯示真實(shí)色彩信息,以圖11和圖12示。對比分析兩種點(diǎn)云濾波結(jié)果,結(jié)果表示:①兩種方法均保留山體整體輪廓特征,上述可視化分析地表地物T1區(qū)域的低矮植被PTD濾波算法結(jié)果基本保留原有特征而改進(jìn)PTD濾波算法剔除部分低矮植被,并未完全清除,植被覆蓋區(qū)域點(diǎn)云出現(xiàn)小面積空洞現(xiàn)象；T2區(qū)域兩種濾波算法均把信號塔剔除；T3區(qū)域建筑物被兩種算法過濾,但PTD濾波算法對蔬菜大棚仍有冗余信息,而改進(jìn)PTD相對來說剔除大部分蔬菜大棚。②實(shí)驗(yàn)區(qū)山體按道路劃分左側(cè)K1區(qū)域和右側(cè)K2、K3區(qū)域PTD濾波算法均誤把山體面點(diǎn)識別為非地面點(diǎn)剔除,從而在山體兩側(cè)出現(xiàn)點(diǎn)云空洞的情況。而改進(jìn)PTD濾波算法對K1、K2和K3區(qū)域的山體地面信息都較為完整保留。

圖11 PTD濾波算法結(jié)果

圖12 改進(jìn)PTD濾波算法結(jié)果

5 結(jié) 論

無人機(jī)DIM點(diǎn)云濾波處理不僅在地物分類、地物單體化提取和地形特征分析中起決定性作用,剔除地物點(diǎn)對后續(xù)生成相關(guān)產(chǎn)品也是核心步驟。對于恐龍谷南緣山區(qū)DIM點(diǎn)云,考慮實(shí)驗(yàn)區(qū)位于滇中高原和山體本身與地面就存在高程差,選擇使用PTD濾波算法對山區(qū)進(jìn)行濾波處理并依據(jù)地形特征點(diǎn)對PTD濾波算法進(jìn)行改進(jìn)試驗(yàn)研究。結(jié)果表明:①明顯地物在被剔除的同時,較為完整保留整個實(shí)驗(yàn)區(qū),明顯地物中山體兩側(cè)低矮植被和山腳蔬菜大棚基本未被剔除,山頂信號塔和山腳建筑物全部清除,但山體部分的地面點(diǎn)易被識別為非地面點(diǎn)而在出現(xiàn)山體K1、K2、K3區(qū)域的空洞現(xiàn)象。②對此,針對恐龍谷南緣山區(qū)復(fù)雜地形,提出以脊-谷交匯地形特征點(diǎn)為改進(jìn)PTD濾波算法的種子點(diǎn),在山體兩側(cè)精細(xì)構(gòu)網(wǎng),山體低矮植被部分清除,相對于PTD濾波算法蔬菜大棚大面積被清除。且山體兩側(cè)未出現(xiàn)明顯點(diǎn)云空洞的現(xiàn)象。

DIM點(diǎn)云通過脊-谷交匯地形特征點(diǎn)改進(jìn)PTD算法對山區(qū)點(diǎn)云濾波相對于經(jīng)典PTD算法更為精準(zhǔn)區(qū)分地面點(diǎn)和非地面點(diǎn),為后續(xù)地形地貌提取和分析提供了空間數(shù)據(jù)支持依據(jù)。但對于復(fù)雜地形和地形起伏度較大區(qū)域,對于地形特征點(diǎn)作為種子點(diǎn)的PTD濾波算后續(xù)研究中對提高提取地形特征點(diǎn)的精確性有待進(jìn)一步加強(qiáng)。