基于ANUDEM與山脊抬升的DEM簡(jiǎn)化

高 翔，徐 柱

（西南交通大學(xué) 遙感信息工程系，四川 成都611756）

為得到粗尺度DEM，需要對(duì)精細(xì)尺度DEM進(jìn)行簡(jiǎn)化［1－2］。現(xiàn)階段，常用的地形簡(jiǎn)化方法有兩類：基于濾波的簡(jiǎn)化和基于地形特征的簡(jiǎn)化［3］。濾波方法常使用移動(dòng)窗口對(duì)DEM進(jìn)行過濾運(yùn)算，如表面簡(jiǎn)化［4］和稀疏采樣［5］等。該類方法雖簡(jiǎn)單高效，但易造成地形的過度平滑從而丟失重要的地形特征。基于特征的簡(jiǎn)化方法首先提取地形結(jié)構(gòu)特征，之后利用得到的特征信息重建粗尺度DEM，因此，基于特征的地形簡(jiǎn)化方法能夠更好地保留地形特征［6－7］。

本文利用ANUDEM（Australia National University DEM）算法基于特征點(diǎn)以及水文要素對(duì)DEM實(shí)施地形簡(jiǎn)化［8］。該方法采用迭代有限微分內(nèi)插技術(shù)和地形強(qiáng)化算法，自動(dòng)去除偽洼地［9－11］，同時(shí)該算法可利用等高線、水文線、洼地、邊界線和海岸線等要素對(duì)插值過程進(jìn)行約束從而得到與現(xiàn)實(shí)地表高度吻合的地形表面。然而，實(shí)驗(yàn)表明，基于特征點(diǎn)與水文的ANUDEM插值方法能夠充分保留地表的水文要素，卻無法顧及山體的結(jié)構(gòu)特征，如山脊線要素。作為地形骨架線，山脊線在很多應(yīng)用中有著比水文線更為重要的作用，如山體單元?jiǎng)澐郑?2－13］與冰川編目［14］。因此，在ANUDEM算法基礎(chǔ)上，提出一種山脊抬升的DEM簡(jiǎn)化方法。該方法將特征點(diǎn)、水文線和山脊線同時(shí)加入到地形簡(jiǎn)化過程中，最終生成基于特征點(diǎn)與骨架線約束的地形簡(jiǎn)化表面。實(shí)驗(yàn)部分，對(duì)不同分辨率下的簡(jiǎn)化DEM進(jìn)行山脊抬升并對(duì)比分析抬升前后的山脊線保留情況。實(shí)驗(yàn)表明，新方法在保證地形參數(shù)準(zhǔn)確性的情況下更好地保留了初始地形的骨架信息

1 方 法

為充分保留地形骨架信息，利用特征點(diǎn)與骨架線進(jìn)行地形簡(jiǎn)化。該方法由4部分組成：①特征點(diǎn)、骨架線的提取；②基于特征點(diǎn)與水文線的ANUDEM內(nèi)插；③簡(jiǎn)化地形的山脊抬升；④山脊抬升前后山脊線保留度統(tǒng)計(jì)分析。本文利用Arc GIS 10．2進(jìn)行DEM、等高線和地形透視圖的展示；利用開源軟件GAT（Geospatial Analysis Tools）完成特征點(diǎn)與水文線的提取；利用Arc Tool box中的Topo to Raster工具完成ANUDEM插值過程，山脊線提取與抬升、地形參數(shù)和骨架線保留度的計(jì)算在IDLE（Python GUI）平臺(tái)中編程實(shí)現(xiàn)。

1．1 提取特征點(diǎn)、骨架線

1．1．1 提取特征點(diǎn)

常用的特征點(diǎn)提取方法有：基于線簡(jiǎn)化的特征點(diǎn)提取算法［15］、VIP算法［16］及Z最大容差法。其中，Z最大容差法是最為常用的特征點(diǎn)提取算法。該方法首先提取一系列點(diǎn)并構(gòu)建TIN；接著向TIN中加入高程誤差（簡(jiǎn)化后DEM與初始DEM高程差）最大的點(diǎn)并重構(gòu)TIN；循環(huán)執(zhí)行第二步直到任意點(diǎn)的高程誤差均小于設(shè)定閾值為止［17］。其中，Z容差閾值與簡(jiǎn)化DEM分辨率的對(duì)應(yīng)關(guān)系如表1所示［18］。

表1 容差值與DEM分辨率的對(duì)應(yīng)關(guān)系 m

1．1．2 提取水文線與山脊線

本文采用水文模擬法提取DEM中的水文線，其一般步驟為①對(duì)初始DEM實(shí)施填洼操作②計(jì)算DEM流向；③基于流向計(jì)算匯水累積量；④基于匯水累積量計(jì)算水文線（本文選取500為累積量判別閾值，因?yàn)楫?dāng)閾值為500時(shí)可獲得較為豐富的水系）；⑤基于流向數(shù)據(jù)，將水文轉(zhuǎn)換為線狀要素。之后，利用文獻(xiàn)［14］設(shè)計(jì)的算法提取DEM中的山脊線，并利用文獻(xiàn)［19］對(duì)地形特征線進(jìn)行優(yōu)化處理。

1．2 山脊抬升

完成ANUDEM插值過程后，對(duì)簡(jiǎn)化表面實(shí)施山脊抬升。若僅增加山脊區(qū)域的高程，山脊與其鄰近區(qū)域?qū)⒋嬖谳^大高差。為減小山脊區(qū)域高程的劇烈變化，本文設(shè)計(jì)了基于最短歐式距離的山脊抬升方案，該方法通過計(jì)算柵格至最近山脊的最短距離來確定當(dāng)前柵格的抬升高度。基于抬升方法的DEM山脊修復(fù)公式為

其中：R為當(dāng)前DEM分辨率，D為某柵格至其最近山脊的最小歐式距離，K為抬升率，H為抬升高程值。

圖1為不同K值和不同H值下，地表抬升高程值與最小歐式距離的函數(shù)曲線，可以發(fā)現(xiàn)：①K值相同時(shí)，H值越大，相同距離下的抬升值越大，且隨著距離的縮小，這一規(guī)律更加明顯；②H值相同時(shí)，K值越大，相同距離下的抬升值越小，其變化趨勢(shì)也越明顯。因此，為減小初始DEM與處理后的DEM間的高程差異，H值不應(yīng)過大，K值則不應(yīng)過小。實(shí)驗(yàn)部分將通過研究地形基本參數(shù)的變化規(guī)律確定H與K的取值。

圖1 抬升高程值與距離的函數(shù)關(guān)系

2 統(tǒng)計(jì)分析

本文選用平均坡度與表面粗糙度對(duì)山脊線抬升前后的地形復(fù)雜度進(jìn)行刻畫。平均坡度ˉS與表面粗糙度K計(jì)算如式（2）、式（3）所示。

其中：S為坡度，A為投影面積，A′為表面面積，i為第i個(gè)柵格單元，n為DEM柵格總數(shù)。為描述山脊線抬升前后地表高程的變化規(guī)律，計(jì)算其高程均方根誤差，如式（4）所示。

圖2所示，式中的L為初始DEM提取出的山脊線原長(zhǎng)（黑色線條）。為計(jì)算地形簡(jiǎn)化后山脊線保留長(zhǎng)度，首先為初始線段創(chuàng)建半徑r＝R／2的緩沖區(qū)，其中R為DEM分辨率；接著統(tǒng)計(jì)簡(jiǎn)化后的骨架線落入緩沖區(qū)內(nèi)的線段長(zhǎng)度L′i；最終基于式（5）獲得保留度［18］。

圖2 骨架線保留度

圖2中實(shí)線為原始DEM骨架線，邊長(zhǎng)L；虛線為簡(jiǎn)化后DEM骨架線；粗線為簡(jiǎn)化后骨架線落入原始DEM骨架線緩沖區(qū)內(nèi)的線段，長(zhǎng)度為L(zhǎng)′i。

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

本文選取GDEM 30 m分辨率DEM進(jìn)行實(shí)驗(yàn)分析，DEM大小及基本參數(shù)如表2所示。

表2 DEM地形基本參數(shù)

由表2可知：實(shí)驗(yàn)所選DEM高差、平均坡度以及表面粗糙度都較大，這表明該區(qū)域具有一定的地表復(fù)雜度，因此適合作為實(shí)驗(yàn)區(qū)域進(jìn)行山脊抬升前后的對(duì)比分析。最終，所選DEM提取的骨架線如圖3所示。

圖3 實(shí)驗(yàn)DEM骨架線提取結(jié)果（黑色細(xì)線為等高線，灰色粗線為水文線，黑色粗線為山脊線）

3．1 山脊抬升對(duì)地形的影響

確定式（1）中山脊抬升的基本參數(shù)K和H，研究不同取值下山脊抬升前后簡(jiǎn)化地表的均方根誤差、平均坡度以及表面粗糙度的變化規(guī)律。實(shí)驗(yàn)數(shù)據(jù)如表3和表4所示。其中，表3為H值相同時(shí)，不同K值下的地形參數(shù)取值；表4為K值相同時(shí)，不同H值下的地形參數(shù)取值。

表3 不同K值下的地形參數(shù)取值

表4 同H值下的地形參數(shù)取值

由表中數(shù)據(jù)可知 與均方根相比 平均坡度與表面粗糙度變化極小，因此山脊抬升不會(huì)顯著影響簡(jiǎn)化地形的基本形態(tài)。圖4進(jìn)一步描述了均方根誤差隨山脊抬升參數(shù)的變化規(guī)律。可以發(fā)現(xiàn)：①H值不變時(shí)，隨著K值的增加，高程均方根誤差呈冪函數(shù)遞減趨勢(shì)，當(dāng)K值小于1時(shí)，均方根誤差迅速減小，K值大于1時(shí)，均方根誤差的變化趨于平緩，如圖4（a）所示，因此實(shí)驗(yàn)部分的K值取1；②K值不變時(shí)，隨著H值的增加，均方根誤差呈線性增長(zhǎng)，如圖4（b）所示，因此為在一定的均方根誤差內(nèi)提升山脊區(qū)域DEM高程，實(shí)驗(yàn)部分的H值取10。

完成抬升參數(shù)的確定后，利用ANUDEM算法與山脊抬升法得到簡(jiǎn)化DEM。其中，50 m、250 m、500 m分辨率下的簡(jiǎn)化DEM如圖5所示。

圖4 均方根誤差與抬升參數(shù)的函數(shù)關(guān)系

圖5 基于ANUDEM與山脊抬升法得到的不同簡(jiǎn)化程度下的DEMs

3．2 山脊保留程度研究

3．2．1 等高線對(duì)比

圖6為基于ANUDEM方法與山脊抬升法得到的50 m分辨率DEM的局部等高線，等高距為100 m。可以發(fā)現(xiàn)：①非山體區(qū)域，山脊抬升前后的等高線近似相同；②山體區(qū)域，山脊抬升前后的等高線具有一定的差異，與原始DEM相比，山脊抬升后的DEM更好地保持了山體特征。圖6（a）丟失了高程為4 600 m的山峰，而圖6（b）保留了這一地形要素，同時(shí)，圖6（b）中高程值為4 700 m的等高線與原始等高線的吻合度也高于圖6（a）。

3．2．2 山脊線保留度對(duì)比

為定量描述抬升前后山脊線保留程度，利用式（5）進(jìn)行保留度計(jì)算，結(jié)果如表5所示。

通過表5數(shù)據(jù)，可以發(fā)現(xiàn)：

1）隨著地形簡(jiǎn)化程度的縮小，特征點(diǎn)百分比逐漸増大，山脊線保留度逐漸增加。如當(dāng)分辨率為500 m時(shí)，抬升后的山脊線保留度僅為43%左右，而當(dāng)分辨率減小至50 m時(shí)，抬升后的山脊線保留度達(dá)到了92%。

圖6 山脊抬升前后DEM等高線差異

2）當(dāng)?shù)匦魏?jiǎn)化程度較高時(shí)，特征點(diǎn)百分比較小，簡(jiǎn)化DEM的數(shù)據(jù)量較小，對(duì)DEM實(shí)施抬升處理并不能夠有效增加山脊線保留度，隨著簡(jiǎn)化程度的降低，山脊抬升能夠一定程度地增加山脊線保留程度。如當(dāng)分辨率等于500 m時(shí)，特征點(diǎn)總數(shù)百分比僅為0．38%，兩者的山脊線保留度均為43%左右，無明顯差異。當(dāng)分辨率減小至50 m時(shí)，與抬升前相比，抬升前后的山脊保留度分別為78%與92%，保留度增加了14%左右。

4 結(jié) 論

本文研究了一種基于ANUDEM與山脊抬升的DEM簡(jiǎn)化方法。該方法首先利用最大容差法提取地形特征點(diǎn)，接著利用水文模擬法提取地表水文線以及山脊線，之后通過ANUDEM插值法獲得簡(jiǎn)化DEM，最終對(duì)簡(jiǎn)化DEM實(shí)施山脊抬升從而獲得顧及骨架線的簡(jiǎn)化地表。

為確定山脊抬升方法的基本參數(shù)，計(jì)算抬升后DEM的均方根誤差、平均坡度以及表面粗糙度。數(shù)據(jù)表明，H值不變時(shí)，均方根誤差隨K值的減小呈冪函數(shù)遞減，當(dāng)K值為2時(shí)，均方根誤差的變化趨于穩(wěn)定；K值不變時(shí)，均方根誤差隨H值的增加呈線性增長(zhǎng)。與均方根誤差相比，山脊抬升對(duì)平均坡度與表面粗糙度的影響較小，因此山脊抬升不會(huì)造成地形復(fù)雜度的較大變化。

為進(jìn)一步研究不同簡(jiǎn)化地形下山脊抬升對(duì)山體特征的保留程度，計(jì)算山脊抬升前后的山脊線保留度。實(shí)驗(yàn)表明，在插值生成的簡(jiǎn)化DEM中實(shí)施山脊抬升能夠一定程度的增加山體結(jié)構(gòu)的保留程度，且隨著地形簡(jiǎn)化程度的降低，山脊抬升法的作用也逐漸明顯。

［1］ ZHANG X，DRAKE N A，WAINWRIHGT J，et al．Co mparison of Slope Esti mates fro m Low Resolution DEMs：Scaling Issues and a Fractal Method f or Their Solution［J］．Earth Surface Processes and Landf or ms，1999，24（9）：763－779．

［2］ 周小軍，王光霞，薛志偉，等．基于DEM化簡(jiǎn)的等高線綜合研究［J］．測(cè)繪工程，2014，23（2）：10－14．

［3］ LI Zhilin．Multi－scale Digital Terrain Modelling and Analysis［C］／／Zhou Q．，Lees B．，Tang G．Advances in Digital Terrain Analysis．Berlin：Springer－Verlag，2008：59－83．

［4］ LI Zhilin，ZHU Q．Digital Elevation model［M］．Wuhan：Wuhan University of Surveying and Mapping Press，2003．

［5］ LI Zhilin．Variation of t he Accuracy of Digital Terrain Models with Sampling Interval［J］．Photogrammetric Record，1992，14（79）：113－128．

［6］ GESCH D B．The Effects of DEM Generalization Methods on Derived Hydr ologic feat ures［C］／／Lowell K，Jaton A．Spatial Accuracy Assess ment：Land Information Uncertainty in Natural Resources．Chelsea：Ann Ar bor Press，1999：255－262．

［7］ AI T，LI J．A DEM Generalization by Minor Valley Branch Detection and Grid Filling［J］．ISPRS Journal of Photogrammetry and Remote Sensing，2010，65（2）：198－207．

［8］ HUTCHINSON M F．ANUDEM Version 4．6．2．http：／／cres．a(chǎn)nu．edu．a(chǎn)u／soft ware／anudem．ht ml as of July 1999．

［9］ 楊勤科，Mcvicar，T．R．，李領(lǐng)濤，等．ANUDEM－專業(yè)化數(shù)字高程模型插值算法及其特點(diǎn)［J］．干旱地區(qū)農(nóng)業(yè)研究，2006，24（3）：36－41．

［10］張彩霞，楊勤科，段建軍．一種高質(zhì)量的數(shù)字高程模型（DEM）建立方法－ANUDEM法［J］．中國農(nóng)業(yè)通報(bào)，2005，21（12）：411－415．

［11］丑述仁，姚志宏，曹佳云，等．基于Hutchinson的DEM建立及質(zhì)量評(píng)價(jià)［J］．地理空間信息，2012，10（3）：127－129．

［12］肖飛，張百平，凌峰．基于DEM的地貌實(shí)體單元自動(dòng)提取方法［J］．地理研究，2008，27（2）：459－467．

［13］徐靜，王春，張耀民，等．規(guī)則格網(wǎng)DEM中平直面狀特征地 形 識(shí) 別 與 提 取［J］．測(cè) 繪 科 學(xué)，2014，39（8）：163－166．

［14］郭萬欽，劉時(shí)銀，余蓬春，等．利用流域邊界和坡向差自動(dòng)提取山脊線［J］．測(cè)繪科學(xué)，2011，36（6）：210－213．

［15］DOUGLAS D H，PEUCKER T K．Algorith ms for the Reduction of the Nu mber of Points Required to Represent a Digitized Line or Its Caricature［J］．Canadian Cartographer，1973，10（2）：112－122．

［16］CHEN Z，GUEVARA J．Systematic Selection of Very Important Points（VIP）fro m Digital Terrain Model for Constructing Triangular Irregular Net works［R］．In：Proc．8thInternational Sy mposium on Computer－Assisted Cartography，Auto Carto8，Balti more，MD，1987，29 March－3 April．50－56．

［17］CHANG K T．Introduction to Geographic Infor mation Systems［M］．Singapore：Mc Graw－Hill，2008．

［18］CHEN Y，WILSON J P，ZHU Q，et al．Comparison of Drainage－constrained Methods f or DEM Generalization［J］．Co mputers and Geosciences，2012，48，41－49．

［19］馮長(zhǎng)強(qiáng)，江南，張巍，等．地形特征線提取后續(xù)問題優(yōu)化處理［J］．測(cè)繪工程，2014，23（3）：28－31．