謝中凱，陳振杰，李飛雪

（1.江蘇省地理信息技術(shù)重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210023）

傳統(tǒng)的由人工解譯地形獲取的等高線數(shù)據(jù)是建立DEM的重要數(shù)據(jù)源[1]。矢量等高線生成規(guī)則格網(wǎng)DEM的方法主要有3種[2,3]：①等高線構(gòu)建TIN法，矢量等高線生成不規(guī)則三角網(wǎng)TIN，然后由TIN內(nèi)插生成DEM；②等高線離散化法，矢量等高線離散化為數(shù)據(jù)點(diǎn)，然后采用一般插值方法（如IDW、Kriging插值法等）內(nèi)插生成DEM[4]；③等高線內(nèi)插法，矢量等高線按照預(yù)定方向直接內(nèi)插生成DEM，或首先進(jìn)行柵格化，然后由柵格化后的等高線內(nèi)插生成DEM。

不同的GIS軟件在利用等高線生成DEM時(shí)采用了不同的方法，ArcGIS采用等高線構(gòu)建TIN法生成DEM，該方法是利用等高線生成DEM的常用方法，但如果在TIN生成時(shí)沒有人工加入地面特征點(diǎn)，生成的TIN會(huì)有平三角形，這樣插值生成的DEM會(huì)產(chǎn)生“平山頂”現(xiàn)象和階梯狀地形。地理資源分析支持系統(tǒng)GRASS是開源GIS軟件[5]，采用等高線內(nèi)插法生成DEM，分為矢量柵格化和柵格插值2個(gè)步驟，其中柵格插值采用線性內(nèi)插方法，通過鄰近的等高線進(jìn)行比例內(nèi)插。本文分析了GRASS軟件中等高線生成DEM算法中存在的問題，并對(duì)算法進(jìn)行了改進(jìn)。

1 GRASS中DEM生成算法存在的問題

1）在矢量柵格化過程中，如果多條等高線通過一個(gè)柵格單元（圖1中為2條），該柵格單元的取值取決于最后進(jìn)行柵格化的等高線的高程值（圖中紅色單元為待插值柵格單元，黑色單元為等高線柵格化后的柵格單元（稱為等高線柵格單元））。

圖1 2條等高線過同一柵格單元

2）在柵格插值過程中，沒有采用最近鄰的2條等高線進(jìn)行內(nèi)插，而是通過四方向Flood Fill算法尋找近鄰等高線。四方向Flood Fill算法的特點(diǎn)是可以避免在搜索時(shí)穿越等高線，具體思路為：對(duì)紅色柵格單元的鄰域按照“下、右、上、左”的順序進(jìn)行搜索判斷（如圖2），對(duì)搜索到的非等高線柵格單元重復(fù)此操作，直到找到2個(gè)不同高程值的等高線柵格單元。當(dāng)算法終止時(shí)（如圖3），搜索結(jié)果為具有不同高程值的藍(lán)色等高線柵格單元，但顯然紫色等高線柵格單元較上側(cè)的藍(lán)色柵格單元距離紅色待插值柵格單元更近，因此在插值時(shí)應(yīng)當(dāng)采用紫色等高線柵格單元。

圖2 四方向Flood Fill算法圖

圖3 四方向Flood Fill搜索結(jié)果圖

3）內(nèi)插生成的DEM存在“平山頂”現(xiàn)象。對(duì)于圖4中的普通情況，A處所在的待插值柵格單元可通過Flood Fill算法找到2條不同高程值的等高線，采用式（1）進(jìn)行線性內(nèi)插。但是對(duì)于圖5中的山頂情況，由于Flood Fill算法搜索時(shí)不能穿越等高線，搜索結(jié)束時(shí)只能找到1條等高線，采用式（2）確定A點(diǎn)的高程值，那么被等高線Z1圍成的山頂將具有相同的高程值，生成的DEM存在“平山頂”現(xiàn)象。

圖4 普通情況

圖5 山頂情況

2 算法改進(jìn)及實(shí)現(xiàn)

針對(duì)GRASS算法存在的3個(gè)問題，分別提出如下解決方案：

1）在矢量柵格化過程中，統(tǒng)計(jì)經(jīng)過每個(gè)柵格單元的等高線的條數(shù)，并記錄其高程值，然后取高程平均值作為該柵格單元的取值。

2）在柵格插值過程中，采用鏈表記錄搜索到的等高線，每條等高線的信息作為一個(gè)鏈表結(jié)點(diǎn)，按距離升序排列，內(nèi)插時(shí)取最近的2條等高線。在搜索過程中，當(dāng)找到第2條等高線后，開始用前2條等高線的最大值作為搜索范圍的限制，減少柵格單元搜索的個(gè)數(shù)，加快搜索速度。

3）在柵格插值過程中，對(duì)于等高線圍成的山頂，采取二輪搜索策略，第二輪搜索時(shí)允許跨越一次等高線，使得山頂采用最近鄰的2條等高線進(jìn)行趨勢(shì)內(nèi)插[6]。在進(jìn)行趨勢(shì)內(nèi)插時(shí)，待插值柵格單元采用式（3）計(jì)算，并且應(yīng)設(shè)定一個(gè)較小的正值，采用式（4）計(jì)算，避免山頂?shù)母叨冗^高。

3 算法實(shí)現(xiàn)及評(píng)價(jià)

等高線生成DEM的改進(jìn)算法利用Visual Studio 2010平臺(tái)和GDAL庫，以C++語言實(shí)現(xiàn)，其中矢量柵格化采用Bresenham算法，柵格插值采用改進(jìn)后的GRASS線性內(nèi)插算法。算法測試選取面積約4 km2的山地區(qū)，原始DEM分辨率為5 m×5 m。首先由原始DEM獲取10 m等高線，然后用不同的算法由等高線重新插值生成DEM。算法評(píng)價(jià)通過3組DEM之間的比較，3組DEM分別為：ArcGIS軟件生成的DEM（未加入人工特征點(diǎn)）、GRASS軟件生成的DEM和本文改進(jìn)算法生成的DEM。

在研究區(qū)內(nèi)隨機(jī)生成500個(gè)隨機(jī)采樣點(diǎn)，進(jìn)行DEM的精度評(píng)估[7]。圖6中，橫軸為原始DEM的采樣值，縱軸為不同算法生成的DEM的采樣值，可以看出3種算法得到的擬合直線斜率均接近1，擬合度R2值達(dá)到0.99以上。但是本文提出的算法擬合直線的斜率要高于其他2種，而截距要小于其他2種，說明改進(jìn)后的算法取得了較好的效果，提高了插值精度。

對(duì)插值效果的評(píng)價(jià)，采用等高線回放法[8]。如圖7所示，紅色線為原始等高線，綠色線為不同算法生成的DEM回放得到的等高線。圖7a、b中回放得到的等高線與實(shí)際等高線相比，均在山頂處出現(xiàn)了失真現(xiàn)象，圖7a中還在等高線曲率較大的地方出現(xiàn)了失真現(xiàn)象，而圖7c中回放得到的等高線與實(shí)際等高線相比，誤差較小。圖7a中出現(xiàn)失真是因?yàn)樵谏巾數(shù)忍厥獾匦翁幧闪似饺切?，?dǎo)致生成的DEM具有“平山頂”現(xiàn)象和階梯狀地形。圖7b中失真是因?yàn)樯巾斕幉逯禃r(shí)直接采用最近等高線賦值，導(dǎo)致生成的DEM具有“平山頂”現(xiàn)象。

圖6 不同算法生成的DEM精度評(píng)價(jià)

圖7 等高線回放比較

在特殊地形處的插值效果方面，選取山頂處DEM進(jìn)行比較。圖8中上圖為不同算法生成的DEM，下圖為相應(yīng)的山體陰影模型。圖8a、b中可以看到明顯的“平山頂”現(xiàn)象，而圖8c中這種現(xiàn)象得到了消除，與原始DEM相接近。

圖8 山頂處DEM生成比較

4 結(jié) 語

對(duì)開源GRASS軟件中等高線生成DEM的算法進(jìn)行改進(jìn)后能提高DEM生成的精度，解決了山體陡峭地區(qū)等高線生成DEM存在的“平山頂”問題。由于現(xiàn)實(shí)中地形是多變的，因此等高線的分布也是各異的，對(duì)于多條等高線經(jīng)過同一柵格單元，下一步的研究可以考慮統(tǒng)計(jì)柵格單元內(nèi)不同高程值的等高線的長度，然后以長度為權(quán)重計(jì)算單元取值。

[1]宋敦江, 岳天祥, 杜正平.由等高線建立DEM的YUEHASM方法研究[J].地球信息科學(xué)學(xué)報(bào), 2009,11(3): 325-332

[2]張凱選, 潘夢(mèng)清, 方輝, 等.利用等高線生成 DEM 方法的研究[J].測繪工程, 2007,16(3): 15-18

[3]楊曉云, 唐咸遠(yuǎn), 梁鑫.基于等高線生成DEM的內(nèi)插算法及其精度分析[J].測繪工程, 2006,15(2):37-39

[4]徐瀟, 譚衢霖, 王浩宇, 等.復(fù)雜地貌地形圖等高線內(nèi)插DEM算法的精度分析[J].遙感信息, 2013,28(6):111-115

[5]吳超.地理資源分析支持系統(tǒng)GRASS的研究與應(yīng)用[D].武漢：華中科技大學(xué)，2005

[6]呂建峰, 劉定生，焦偉利, 等.DEM生成算法并行化研究[J].中國圖像圖形學(xué)報(bào),2002, 7(5):506-512

[7]寇程, 柯長青.地形平坦地區(qū)DEM生成算法的比較研究[J].測繪與空間地理信息,2013, 36(7): 33-37

[8]唐新明, 林宗堅(jiān), 吳嵐.基于等高線和高程點(diǎn)建立DEM的精度評(píng)價(jià)方法探討[J].遙感信息,1999,14(3):7-10