恐龍谷南緣山區DIM點云的PTD濾波改進試驗研究

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

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

1 引 言

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

對于祿豐恐龍谷南緣山區特殊地形地貌,采用低空無人機測量技術獲取的影像數據,并構建DIM點云,首先考慮山頂、山體兩側和地面存在高程差,選擇漸進加密三角網濾波算法進行處理。其次針對實驗區特殊地理位置和地貌形態,改進傳統漸進加密三角網濾波算法,最后對兩種濾波算法處理結果進行直觀分析。較為完整得到實驗區地面點數據可為后續解譯地貌特征信息和提取地形特征信息奠定基礎。

2 PTD改進方法研究

2.1 DIM點云PTD濾波方法

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

圖1 判斷地面點

2.2 PTD算法的改進技術

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

為此需要不斷剔除最遠點和增加點構建移動曲面,確保點云規則傳遞并覆蓋整個實驗區,如3圖示,主要進行:①把實驗區固定格網按建筑物最大尺寸劃分；②確定實驗區最外圍格網高程最低點,以此為地面種子點；③點A、B、C、D是研究區四個角格網中心,以此構網,域內全部點為中心建立移動格網；④如圖4示,進行多尺度迭代把移動格網內最低高程點作為地面種子點；⑤剔除重復點,并和格網外地面種子點建立實驗區地面種子點。

圖2 表示待定點識別

圖3 移動格網構建

圖4 格網多尺度迭代

為避免格網缺失產生的錯誤分類,格網ABCD外部不再構建移動格網。每個移動格網均依照實驗區建筑物最大尺寸劃分,既保證了所有格網均有地面點,又確保獲取地面種子點無需進一步優化,直接進入三角網迭代判斷的步驟。

3 試驗區及試驗數據采集

3.1 試驗區

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

圖5 實驗區概況

3.2 試驗數據采集

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

4 實驗分析

4.1 基于DIM點云可視化分析

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

圖6 數字表面模型(DSM)

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

圖7 正射影像圖(DOM)

4.2 基于PTD改進的點云濾波處理成效對比分析

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

實驗通過“試錯法”調整閾值,選定最大地形坡度角為88°,迭代角度閾值設置為6°,迭代距離設置為1.4 m。70個初試種子點選取如圖8(a)示,種子點除了邊界點,其余均勻分布于實驗區,由70個種子點構建TIN如圖8(b)示,構建三角網較為規整。

(a) 濾波種子點

(b) 種子點構建不規則三角網

實驗區采集DIM點云有17271784個點,由初試種子點構建TIN進行PTD濾波,剔除非地面點并進行冗余數據抽希,最后保留地面點5723274個點。

4.2.2 改進PTD濾波處理技術應用

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

圖9 主要脊-谷線示意

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

(1)

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

表1 地形特征點提取精度評價

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

圖10 (a)地形特征點

圖10 (b)特征點構建不規則三角網

4.2.3 改進PTD去噪結果對比分析

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

圖11 PTD濾波算法結果

圖12 改進PTD濾波算法結果

5 結 論

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

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