基于ANUDEM與山脊抬升的DEM簡化

高 翔，徐 柱

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

為得到粗尺度DEM，需要對精細(xì)尺度DEM進(jìn)行簡化［1－2］。現(xiàn)階段，常用的地形簡化方法有兩類：基于濾波的簡化和基于地形特征的簡化［3］。濾波方法常使用移動窗口對DEM進(jìn)行過濾運算，如表面簡化［4］和稀疏采樣［5］等。該類方法雖簡單高效，但易造成地形的過度平滑從而丟失重要的地形特征?；谔卣鞯暮喕椒ㄊ紫忍崛〉匦谓Y(jié)構(gòu)特征，之后利用得到的特征信息重建粗尺度DEM，因此，基于特征的地形簡化方法能夠更好地保留地形特征［6－7］。

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

1 方 法

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

1．1 提取特征點、骨架線

1．1．1 提取特征點

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

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

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

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

1．2 山脊抬升

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

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

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

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

2 統(tǒng)計分析

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

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

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

圖2 骨架線保留度

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

3 實驗與分析

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

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

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

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

3．1 山脊抬升對地形的影響

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

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

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

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

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

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

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

3．2 山脊保留程度研究

3．2．1 等高線對比

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

3．2．2 山脊線保留度對比

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

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

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

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

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

4 結(jié) 論

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

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

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

［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化簡的等高線綜合研究［J］．測繪工程，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ù)字高程模型插值算法及其特點［J］．干旱地區(qū)農(nóng)業(yè)研究，2006，24（3）：36－41．

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

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

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

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

［14］郭萬欽，劉時銀，余蓬春，等．利用流域邊界和坡向差自動提取山脊線［J］．測繪科學(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］馮長強(qiáng)，江南，張巍，等．地形特征線提取后續(xù)問題優(yōu)化處理［J］．測繪工程，2014，23（3）：28－31．