謝中凱,陳振杰,李飛雪
(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í)沒(méi)有人工加入地面特征點(diǎn),生成的TIN會(huì)有平三角形,這樣插值生成的DEM會(huì)產(chǎn)生“平山頂”現(xiàn)象和階梯狀地形。地理資源分析支持系統(tǒng)GRASS是開(kāi)源GIS軟件[5],采用等高線內(nèi)插法生成DEM,分為矢量柵格化和柵格插值2個(gè)步驟,其中柵格插值采用線性內(nèi)插方法,通過(guò)鄰近的等高線進(jìn)行比例內(nèi)插。本文分析了GRASS軟件中等高線生成DEM算法中存在的問(wèn)題,并對(duì)算法進(jìn)行了改進(jìn)。
1)在矢量柵格化過(guò)程中,如果多條等高線通過(guò)一個(gè)柵格單元(圖1中為2條),該柵格單元的取值取決于最后進(jìn)行柵格化的等高線的高程值(圖中紅色單元為待插值柵格單元,黑色單元為等高線柵格化后的柵格單元(稱為等高線柵格單元))。
圖1 2條等高線過(guò)同一柵格單元
2)在柵格插值過(guò)程中,沒(méi)有采用最近鄰的2條等高線進(jìn)行內(nèi)插,而是通過(guò)四方向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處所在的待插值柵格單元可通過(guò)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 山頂情況
針對(duì)GRASS算法存在的3個(gè)問(wèn)題,分別提出如下解決方案:
1)在矢量柵格化過(guò)程中,統(tǒng)計(jì)經(jīng)過(guò)每個(gè)柵格單元的等高線的條數(shù),并記錄其高程值,然后取高程平均值作為該柵格單元的取值。
2)在柵格插值過(guò)程中,采用鏈表記錄搜索到的等高線,每條等高線的信息作為一個(gè)鏈表結(jié)點(diǎn),按距離升序排列,內(nèi)插時(shí)取最近的2條等高線。在搜索過(guò)程中,當(dāng)找到第2條等高線后,開(kāi)始用前2條等高線的最大值作為搜索范圍的限制,減少柵格單元搜索的個(gè)數(shù),加快搜索速度。
3)在柵格插值過(guò)程中,對(duì)于等高線圍成的山頂,采取二輪搜索策略,第二輪搜索時(shí)允許跨越一次等高線,使得山頂采用最近鄰的2條等高線進(jìn)行趨勢(shì)內(nèi)插[6]。在進(jìn)行趨勢(shì)內(nèi)插時(shí),待插值柵格單元采用式(3)計(jì)算,并且應(yīng)設(shè)定一個(gè)較小的正值,采用式(4)計(jì)算,避免山頂?shù)母叨冗^(guò)高。
等高線生成DEM的改進(jìn)算法利用Visual Studio 2010平臺(tái)和GDAL庫(kù),以C++語(yǔ)言實(shí)現(xiàn),其中矢量柵格化采用Bresenham算法,柵格插值采用改進(jìn)后的GRASS線性內(nèi)插算法。算法測(cè)試選取面積約4 km2的山地區(qū),原始DEM分辨率為5 m×5 m。首先由原始DEM獲取10 m等高線,然后用不同的算法由等高線重新插值生成DEM。算法評(píng)價(jià)通過(guò)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種,說(shuō)明改進(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生成比較
對(duì)開(kāi)源GRASS軟件中等高線生成DEM的算法進(jìn)行改進(jìn)后能提高DEM生成的精度,解決了山體陡峭地區(qū)等高線生成DEM存在的“平山頂”問(wèn)題。由于現(xiàn)實(shí)中地形是多變的,因此等高線的分布也是各異的,對(duì)于多條等高線經(jīng)過(guò)同一柵格單元,下一步的研究可以考慮統(tǒng)計(jì)柵格單元內(nèi)不同高程值的等高線的長(zhǎng)度,然后以長(zhǎng)度為權(quán)重計(jì)算單元取值。
[1]宋敦江, 岳天祥, 杜正平.由等高線建立DEM的YUEHASM方法研究[J].地球信息科學(xué)學(xué)報(bào), 2009,11(3): 325-332
[2]張凱選, 潘夢(mèng)清, 方輝, 等.利用等高線生成 DEM 方法的研究[J].測(cè)繪工程, 2007,16(3): 15-18
[3]楊曉云, 唐咸遠(yuǎn), 梁鑫.基于等高線生成DEM的內(nèi)插算法及其精度分析[J].測(cè)繪工程, 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].中國(guó)圖像圖形學(xué)報(bào),2002, 7(5):506-512
[7]寇程, 柯長(zhǎng)青.地形平坦地區(qū)DEM生成算法的比較研究[J].測(cè)繪與空間地理信息,2013, 36(7): 33-37
[8]唐新明, 林宗堅(jiān), 吳嵐.基于等高線和高程點(diǎn)建立DEM的精度評(píng)價(jià)方法探討[J].遙感信息,1999,14(3):7-10