一種顧及地形特征的山區(qū)LiDAR高精度DEM提取算法

彭檢貴 邢元軍 宋亞斌 王宗躍

摘要：針對目前常規(guī)DEM插值算法無法很好地表達斷裂區(qū)域地形細節(jié)的問題，提出了一種基于地形特征約束的高精度DEM插值算法。首先利用離散點構建無約束TIN，然后利用邊緣檢測算法提取地形斷裂特征，并嵌入地形斷裂線作為約束，構建約束TIN，在此基礎上制作高精度DEM。最后，利用黃土高坡數據進行了DEM提取實驗。結果證明，該算法較傳統算法在地形細節(jié)的保留上更具優(yōu)勢。

關鍵詞：地形特征；數字高程模型；激光雷達；約束三角網

DOIDOI：10.11907/rjdk.171700

中圖分類號：TP312文獻標識碼：A文章編號：16727800（2017）010002704

0引言

數字高程模型（Digital Elevation Model，簡稱DEM）是用一組有序數值陣列形式表示地面高程的一種實體地面模型[1]，是數字地形模型（Digital Terrain Model，簡稱DTM）的一個分支[2]。DEM是描述包括高程在內的各種地貌因子的基礎，在測繪、水文、氣象、地貌、地質、土壤、工程建設、通訊、軍事等國民經濟和國防建設領域有著廣泛應用[3]。DEM數據源的獲取方式多種多樣，如：傳統地面測量、攝影測量、干涉合成孔徑雷達（InSAR）和機載激光雷達（Light Detection and Ranging，簡稱LiDAR）等[4]。地面測量精度高，但耗時又費力，不適合大規(guī)模應用；被動方式的攝影測量在植被覆蓋區(qū)域和不利氣象條件下效果不理想[5]；機載激光雷達和InSAR技術由于具有較高精度，并能進行較低成本的大范圍應用，使其成為了工程應用的首選[6]。

機載激光雷達集激光掃描系統、全球定位系統（GPS）和慣性導航系統（INS）三種技術于一身，無需大量地面控制點，即能快速準確地獲取地表高密度、高精度的高程信息[7]，已逐漸成為制作高精度DEM產品的一種優(yōu)質數據源。從機載LiDAR點云數據到DEM產品，需要一定處理流程，其中最重要的步驟是點云濾波和DEM插值。濾波是指從點云中分離地面點、非地面點（包括人工目標（如建筑物、橋梁）和自然地物表面（如：灌木、樹木等）的過程。作為制作DEM的重要步驟，濾波一直是國內外學者的研究重點。目前，已形成了多種代表性算法，如數學形態(tài)學的濾波方法[810]、基于分割后拓撲重建的濾波方法[11，12]、基于坡度的濾波算法[1314]、漸近三角網加密算法（PTD）[1517]等。通過ISPRS第三工作組對幾種常用濾波算法的實驗與分析[18]，Axelsson認為漸近三角網加密算法在地形保留和誤差抑制上做的最好。而專門針對DEM插值的研究并不多，目前常用的插值算法有移動曲面擬合法、徑向基函數法、反距離加權法、克里金插值法、線性三角網插值算法等[19]。

以上幾種算法是在地形連續(xù)的假設上利用數學公式擬合地形，在點云密度較高的地區(qū)能取得較好效果，但在地形起伏較大、點云分布稀疏的區(qū)域則存在較多的信息損失。在植被密布的山區(qū)，LiDAR穿透力有限，經過濾波處理后，植被下方的地面點云密度可能非常低，此時傳統插值方法已不再適用，需要顧及地形的骨架信息，盡可能地保留地形細節(jié)。

針對上述問題，本文提出一種顧及地形骨架的山區(qū)LiDAR高精度DEM提取方法。首先，利用普適性較好的漸近加密三角網算法（PTD）濾波獲取全部地面腳點，同時獲取無約束TIN；然后利用邊緣檢測算法提取地形骨架線（包括山脊線、山谷線以及斷裂線），將地形骨架線作為約束條件，構建約束TIN，并基于約束TIN制作高精度DEM；最后，利用黃土高坡的點云進行了DEM提取實驗。實驗結果證明，本文提出的插值算法在地形起伏較大的山區(qū)能夠獲得精度更高、更加貼近真實地形的DEM。

1顧及地形特征的山區(qū)LiDAR高精度DEM提取

根據約束條件的有無，三角剖分可分成常規(guī)三角剖分（Delaunay Triangulation，簡稱DTIN）和約束三角剖分（Constrained delaunay Triangulation，簡稱CDTIN）。利用離散點云構建三角網時，不僅對三角形的形狀有要求，而且離散數據本身也會影響三角網的局部合理性。通過對一些特殊地物或地形的點（如斷裂、山脊、山谷、堤壩等）進行組合，對TIN進行強制性約束，從而使構建的三角網更符合實際地形。因此，基于約束三角網內插的DEM也更加貼近真實地形。本文在無約束TIN中嵌入地形結構線，構建約束TIN，最后基于約束TIN制作高精度DEM。

1.1基于PTD的地面點提取與非約束三角網構建

Axelsson提出的PTD算法以三角網為基礎，通過由粗到細的過程，利用一定的光滑條件，逐步獲取精密的地面點。由于該算法具有較好的普適性，已在商業(yè)化的LiDAR數據處理軟件（如芬蘭的Terra Solid、中國的LiDAR_Suite）中實現，算法的基本步驟如下：①對LiDAR點云構建格網索引，所需的參數一般需要人為設定；②對于格網的每個分塊，搜索其最低點作為初始地面點，并構建稀疏的地形TIN；③根據待判斷點與所在三角形的角度、距離、鏡像點原則等條件作判斷（見圖1），確定余下LiDAR腳點是否能夠加入三角網，若滿足條件，則將其加入三角網中；④重復步驟③，直到所有點都已被判定為地面點或非地面點。

圖1PTD算法加密過程

PTD作為經典的濾波算法之一，已得到了廣泛研究，其具體細節(jié)在文獻[15][17]中有著非常翔實的描述，本文不再贅述。本文選擇PTD算法濾波獲取地面腳點，首先考慮其極好的普適性，在地形起伏大的山區(qū)也能有效保留地形骨架并較好地抑制噪聲；另一方面，PTD算法在加密濾波的同時也生成三角網，極大簡化了后續(xù)插值工作。

1.2嵌入地形特征構建約束TIN

地形特征線，主要包括山脊線、山谷線以及斷裂線。作為地形的骨架，地形特征線不僅決定了地形地貌的幾何形態(tài)和基本走勢，也有特定的物理意義。本文利用國產激光雷達數據處理軟件LiDAR_Suite基于離散的地面點云提取斷裂線、山脊線和山谷線，將上述3種地形特征線作為約束條件嵌入TIN模型，構建約束Delaunay三角網。endprint

約束邊嵌入的算法很多，如分割合并法[20]、Shell三角化法[21]以及對角線交換算法[22]等。綜合考慮實驗難度以及插入效率，本文選擇對角線交換算法構建約束TIN。對角線交換法的基本思想是，從起始點出發(fā)，判斷每一條對角線的可交換性，若可以則交換，反之則繼續(xù)迭代，直到約束邊完整嵌入三角網為止。為了避免在影響域比較復雜時，判斷和交換也隨之變得復雜從而導致迭代失敗的情況，本文采取m+2多邊形對角線交換法，通過約束邊將影響域分成2個區(qū)域，然后刪除這兩個多邊形區(qū)域內的三角形，重新在該區(qū)域生成三角形，具體流程如下：

（1）假設有顏色區(qū)域為斷裂區(qū)域（見圖2），BL（淺灰色部分）表示高程較高的平坦區(qū)域，BR（深灰色部分）表示斷裂線下方區(qū)域，L1L2為斷裂線段（約束線段）。將約束線段的各個頂點看成離散點，將約束線段頂點加入到獲得的CDT中。利用LOP算法對三角網進行調整，使之滿足LiLj（i，j

（2）根據約束邊確定影響域L。假設正在處理的約束線段為L1L2，則與L1L2相交的三角形構成的區(qū)域稱為約束線段L1L2的影響域MT={T＼-1，T＼-2，…，T＼-n}，其中MT中三角形的外邊組成了該影響域的影響多邊形B={L1，a，b，c，……，L2}（圖2中有顏色的區(qū)域）。該影響多邊形有3個特點：①B是簡單多邊形，L1L2為B的一條對角線，將B分成BL和BR兩部分，并且BL和BR也必須是簡單多邊形（見圖2）；②BL和BR也可以進行Delaunay三角剖分；③ 對于任意Lk，假如Lk是距離L1L2最近的點（Lk≠L1，Lk≠L2），可確定LiLk∈B，LkLj∈B。刪除L1L2上下兩側的多邊形BL和BR內的三角形，重新進行Delaunay 三角剖分，從而將約束線段L1L2嵌入到三角網中。以BL為例，實現步驟如下：①建立一個空的堆棧stack，首先將L1L2的邊放入stack中；②在stack中彈出一邊，設為LcLiLj，遍歷BL頂點，尋找離LiLj最近的點Lc（Lc≠Li，Lc≠Lj），形成三角形LcLiLj；③假如LcLi（或LcLj）是BL的一條邊，則不入堆棧；反之則將LcLi（或LcLj）放入堆棧stack中；④若堆棧不為空，則返回步驟②，否則針對BL的操作結束；⑤ 使用LOP算法對BL和BR的剖分三角形進行局部優(yōu)化處理，使之滿足Delaunay三角形的兩個基本性質。

如圖2（a）所示，假如直接基于無約束三角網進行線性插值（具體細節(jié)見1.3節(jié)），由于三角形穿越斷裂線，BL部分插值得到的高程將明顯低于真實高程，而BR部分的高程則可能大于真實值。而在嵌入地形特征線構建約束三角網之后，地形結構線將不再穿越三角網（見圖2（c）），線性插值獲取的高程值更加貼近真實地形。

1.3基于約束TIN提取顧及地形特征的DEM

建立Delaunay三角網（TIN）之后，即可基于TIN計算區(qū)域內任意一點的高程，在點密度均勻的情況下，可以通過線性內插的方式獲取各網點高程，即令三角形三點確定的傾斜平面作為該微小面元的地表。算法原理如圖3所示。

對于待插值的格網點P（x，y），首先根據P點的平面坐標確定P隸屬的分塊格網號i，然后依次計算待插值點P與格網內其它點距離的平方：

D2i=（x-xi）2+（y-yi）2（2）

求取距離最小的點，設為Q1，然后依次取出Q1為頂點的三角形，判斷點P是否位于該三角形內。若是，停止判斷，反之則繼續(xù)判斷。若Q1為頂點的三角形都不包含點P，則繼續(xù)判斷距離次近的點Q2，如此迭代直到找到點P的外包三角形。

第三步：若P（x，y）的外包三角形為ΔQ1Q2Q3，三頂點的坐標分別是（x1，y1，z1）、（x2，y2，z2）和（x3，y3，z3），然后通過線性插值法確定P點的高程：

Z=Z1-（x-x1）（y21z31-y31z21）+（y-y1）（z21x31-z31x21）x21y31-x31y21（3）

2實驗與結果分析

為了考察本算法的實際運行效果，本文基于國產激光雷達平臺LiDAR_Suite進行二次開發(fā)，在VC++的環(huán)境下，實現了本文算法，并進行了實驗與結果分析。實驗數據選自黃土高坡地區(qū)，該區(qū)域地形復雜、地表形態(tài)豐富。數據采集于2011年9月，點云總數是1 562 831，點云平均間距約0.7m，垂直誤差約0.15m。

分析圖6可知，基于無約束TIN線性插值獲取的DEM在地形平緩的區(qū)域，能夠較好地保留地表細節(jié)，但在斷裂區(qū)域附近，存在較明顯的“平滑”現象。與此同時，如圖7所示，基于約束TIN插值獲取的DEM，其地表斷裂處的地形細節(jié)都保留得更為完整。從視覺上分析，本文插值算法比傳統不顧及地形骨架的方法更適用于山區(qū)DEM的插值。

為了更好地評價本文插值算法的效果，將基于GPS RTK量測獲得的斷裂地形附近的三維坐標作為檢查點，定量評價基于移動曲面擬合、反距離加權法、線性三角網（無約束）、線性三角網（有約束）獲取的DEM精度，其結果如表1所示。

表1中，A＼-1～A＼-6為地形較平緩區(qū)域所量測的檢查點，A＼-7～A＼-11為斷裂地形附近量測的檢查點。分析表1可知，幾種傳統插值算法獲取的DEM精度差異并不大，中誤差在0.5～0.6m之間，約為點云高程誤差的3～4倍，滿足國家機載LiDAR數據后處理規(guī)范山區(qū)1：2 000 DEM的成圖要求。另一方面，地形突變區(qū)域（A＼-7～A＼-11）的誤差明顯大于平坦地區(qū)（A＼-1～A＼-6）的誤差，可見傳統插值算法在地形細節(jié)保真方面有較大的提高空間。相對而言，基于約束TIN插值獲取的DEM，在地形平緩區(qū)域的誤差和傳統插值結果接近，在地形突變結構附近區(qū)域的誤差卻明顯小于傳統插值結果，可見顧及地形骨架信息之后，DEM的提取精度有了明顯提升。

3結語

本文針對傳統DEM插值算法無法很好地表達地形起伏大區(qū)域地表細節(jié)的問題，提出了一種顧及地形骨架的高精度DEM插值算法。首先利用離散點構建無約束TIN，然后利用邊緣檢測算法提取地形骨架線，并嵌入地形斷裂線作為約束，構建約束TIN，在此基礎上利用約束TIN制作高精度DEM。實驗結果證明，本文算法提取的山區(qū)DEM較傳統算法精度更高，地形細節(jié)的保留更為逼真。

參考文獻參考文獻：

[1]湯國安，李發(fā)源，劉學軍.數字高程模型教程[M].北京：科學出版社，2010.

[2]周啟鳴，劉學軍.數字地形分析[M].北京：科學出版社，2006.

[3]張瑞軍，楊武年，劉漢湖，等.數字高程模型（DEM）的構建及其應用[J].工程勘察，2005（5）：6164.

[4]MAGUYA A S， JUNTTILA V， KAURANNE T. Adaptive algorithm for large scale dtm interpolation from lidar data for forestry applications in steep forested terrain[J]. Isprs Journal of Photogrammetry & Remote Sensing， 2013，85（11）：7483.

[5]高廣，馬洪超，張良，等.顧及地形斷裂線的LiDAR點云濾波方法研究[J].武漢大學學報：信息科學版，2015，40（4）：474478.

[6]WILSON J P. Digital terrain modeling[J]. Geomorphology， 2012， 137（1）：107121.

[7]黃先鋒，李卉，王瀟，等.機載LiDAR數據濾波方法評述[J].測繪學報，2009，38（5）：466460.

[8]ZHANG K， CHEN SC， WHITMAN D，et al. A progressive morphological filter for removing nonground measurements from airborne LIDAR data[J]. Geoscience and Remote Sensing， IEEE Transactions on， 2003，41（4）：872882.

[9]CHEN Q， GONG P， BALDOCCHI D， et al. Filtering airborne laser scanning data with morphological methods[J]. Photogrammetric Engineering and Remote Sensing， 2007，73（2）：175185.

[10]PINGEL TJ， CLARKE KC， MCBRIDE WA. An improved simple morphological filter for the terrain classification of airborne LIDAR data[J]. ISPRS Journal of Photogrammetry and remote Sensing， 2013，77（1）：2130.

[11]SITHOLE G， VOOR GEODESIE NC. Segmentation and classification of airborne laser scanner data[M]. Nederlandse Commissie voor Geodesie， 2005.

[12]彭檢貴，馬洪超，高廣，等.利用機載LiDAR點云數據提取城區(qū)道路[J].測繪通報，2012（9）：1619.

[13]SITHOLE G. Filtering of laser altimetry data using a slope adaptive filter[J]. International Archives of Photogrammetry Remote Sensing and Spatial Information Sciences， 2011，34（3/W4）：203210.

[14]SUSAKI J. Adaptive slope filtering of airborne LiDAR data in urban areas for digital terrain model （DTM） generation[J]. Remote Sensing， 2012，4（4）：18041819.

[15]AXELSSON P. DEM generation from laser scanner data using adaptive TIN models[J].International Archives of Photogrammetry and Remote Sensing， 2000，33（B4/1）：111118.

[16]左志權，張祖勛，張劍清.知識引導下的城區(qū)LiDAR點云高精度三角網漸進濾波方法[J].測繪學報，2012，42（2）：246251.

[17]隋立春，張熠斌，張碩，等.基于漸進三角網的機載LiDAR點云數據濾波[J].武漢大學學報：信息科學版，2011，36（10）：11591163.

[18]SITHOLE G， VOSSELMAN G. Experimental comparison of filter algorithms for bareearth extraction from airborne laser scanning point clouds[J]. ISPRS Journal of Photogrammetry and remote Sensing， 2004，59（12）：85101.

[19]張良，馬洪超，鄔建偉.聯合機載LiDAR數據和潮汐數據自動提取潮位線[J].遙感學報，2012，16（2）：405416.

[20]LEE D T， B J SCHACHTER. Two algorithms for constructing a delaunay triangulation[J]. International Journal of Parallel Programming， 1980，9（3）：219242.

[21]PIEGL L A， RICHARD A M. Algorithm and data structure for triangulating multiply connected polygonal domains[J]. Computers & Graphics， 1993，17（93）：563574.

[22]劉學軍，龔健雅.約束數據域的Delaunay三角剖分與修改算法[J].測繪學報，2001，30（1）：8288.

責任編輯（責任編輯：黃健）endprint