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

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

關(guān)鍵詞：地形特征；數(shù)字高程模型；激光雷達(dá)；約束三角網(wǎng)

DOIDOI：10.11907/rjdk.171700

中圖分類號(hào)：TP312文獻(xiàn)標(biāo)識(shí)碼：A文章編號(hào)：16727800（2017）010002704

0引言

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

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

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

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

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

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

1.1基于PTD的地面點(diǎn)提取與非約束三角網(wǎng)構(gòu)建

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

圖1PTD算法加密過(guò)程

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

1.2嵌入地形特征構(gòu)建約束TIN

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

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

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

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

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

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

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

對(duì)于待插值的格網(wǎng)點(diǎn)P（x，y），首先根據(jù)P點(diǎn)的平面坐標(biāo)確定P隸屬的分塊格網(wǎng)號(hào)i，然后依次計(jì)算待插值點(diǎn)P與格網(wǎng)內(nèi)其它點(diǎn)距離的平方：

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

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

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

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

2實(shí)驗(yàn)與結(jié)果分析

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

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

為了更好地評(píng)價(jià)本文插值算法的效果，將基于GPS RTK量測(cè)獲得的斷裂地形附近的三維坐標(biāo)作為檢查點(diǎn)，定量評(píng)價(jià)基于移動(dòng)曲面擬合、反距離加權(quán)法、線性三角網(wǎng)（無(wú)約束）、線性三角網(wǎng)（有約束）獲取的DEM精度，其結(jié)果如表1所示。

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

3結(jié)語(yǔ)

本文針對(duì)傳統(tǒng)DEM插值算法無(wú)法很好地表達(dá)地形起伏大區(qū)域地表細(xì)節(jié)的問(wèn)題，提出了一種顧及地形骨架的高精度DEM插值算法。首先利用離散點(diǎn)構(gòu)建無(wú)約束TIN，然后利用邊緣檢測(cè)算法提取地形骨架線，并嵌入地形斷裂線作為約束，構(gòu)建約束TIN，在此基礎(chǔ)上利用約束TIN制作高精度DEM。實(shí)驗(yàn)結(jié)果證明，本文算法提取的山區(qū)DEM較傳統(tǒng)算法精度更高，地形細(xì)節(jié)的保留更為逼真。

參考文獻(xiàn)參考文獻(xiàn)：

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

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

[3]張瑞軍，楊武年，劉漢湖，等.數(shù)字高程模型（DEM）的構(gòu)建及其應(yīng)用[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點(diǎn)云濾波方法研究[J].武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2015，40（4）：474478.

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

[7]黃先鋒，李卉，王瀟，等.機(jī)載LiDAR數(shù)據(jù)濾波方法評(píng)述[J].測(cè)繪學(xué)報(bào)，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]彭檢貴，馬洪超，高廣，等.利用機(jī)載LiDAR點(diǎn)云數(shù)據(jù)提取城區(qū)道路[J].測(cè)繪通報(bào)，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]左志權(quán)，張祖勛，張劍清.知識(shí)引導(dǎo)下的城區(qū)LiDAR點(diǎn)云高精度三角網(wǎng)漸進(jìn)濾波方法[J].測(cè)繪學(xué)報(bào)，2012，42（2）：246251.

[17]隋立春，張熠斌，張碩，等.基于漸進(jìn)三角網(wǎng)的機(jī)載LiDAR點(diǎn)云數(shù)據(jù)濾波[J].武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，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]張良，馬洪超，鄔建偉.聯(lián)合機(jī)載LiDAR數(shù)據(jù)和潮汐數(shù)據(jù)自動(dòng)提取潮位線[J].遙感學(xué)報(bào)，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]劉學(xué)軍，龔健雅.約束數(shù)據(jù)域的Delaunay三角剖分與修改算法[J].測(cè)繪學(xué)報(bào)，2001，30（1）：8288.

責(zé)任編輯（責(zé)任編輯：黃?。〆ndprint