孫嘉駿,黃天勇,陳 科

（1.中國水電顧問集團(tuán) 昆明勘測(cè)設(shè)計(jì)研究院有限公司 測(cè)繪地理信息分院,云南 昆明 650041）

DEM是通過有限的地形高程數(shù)據(jù)實(shí)現(xiàn)對(duì)地形曲面的數(shù)字化模擬或是地形表面形態(tài)的數(shù)字化表示的量化模型[1]。近年來，隨著空間數(shù)據(jù)基礎(chǔ)設(shè)施的建設(shè)和“數(shù)字地球”戰(zhàn)略的實(shí)施，加快了DEM與地理信息系統(tǒng)、遙感等的一體化進(jìn)程，為DEM的應(yīng)用開辟了更廣闊的天地。

1 DEM概述

1.1 DEM的含義

DEM表示三維向量有限序列，是對(duì)地球表面地形地貌的一種離散的數(shù)字表達(dá)，也是測(cè)繪學(xué)中最常用的地形表示模型，用函數(shù)形式描述為：iiii

式中，（Xi,Yi）是平面坐標(biāo)；Zi是（Xi,Yi）對(duì)應(yīng)的高程。當(dāng)該序列中各平面向量的平面位置呈規(guī)則格網(wǎng)排列時(shí)，其平面坐標(biāo)可省略，此時(shí)DEM就簡化為一維向量序列{Zi, i=1,2,3,…,n}[1]。

1.2 DEM的表示模型及特點(diǎn)

隨著對(duì)DEM需求的不斷增加，對(duì)DEM數(shù)據(jù)的多樣性、實(shí)用性及精度提出了更高的要求。同時(shí)由于其數(shù)據(jù)來源不同，導(dǎo)致數(shù)據(jù)存儲(chǔ)結(jié)構(gòu)的差別。DEM數(shù)據(jù)結(jié)構(gòu)表示模型主要包括：

1）不規(guī)則離散點(diǎn)型。由不規(guī)則離散點(diǎn)表示的DEM是將連續(xù)地球表面形態(tài)離散成在某一個(gè)區(qū)域D上的以Xi、Yi、Zi三維坐標(biāo)形式存儲(chǔ)的高程點(diǎn)Zi((Xi,Yi)ID)的集合[2]。離散點(diǎn)DEM往往是通過測(cè)量直接獲取的地球表面的原始特征點(diǎn)數(shù)據(jù)。這些特征點(diǎn)之間相互獨(dú)立，彼此間沒有任何聯(lián)系，它是DEM中最簡單的數(shù)據(jù)組織形式。

2）規(guī)則格網(wǎng)型。規(guī)則格網(wǎng)DEM是指在水平方向和垂直方向按等間隔方式記錄格網(wǎng)點(diǎn)的坐標(biāo)，格網(wǎng)點(diǎn)的高程Z表示地形，格網(wǎng)點(diǎn)的平面坐標(biāo)隱含在行列號(hào)中，故適宜用矩陣形式進(jìn)行存儲(chǔ)，即按行(或列)逐一記錄每一個(gè)格網(wǎng)單元的高程值。通常，為了壓縮柵格DEM的冗余數(shù)據(jù)，可采取用游程編碼或四叉樹編碼方法對(duì)數(shù)據(jù)進(jìn)行處理；為準(zhǔn)確表示地形細(xì)部，可采用附加地形特征數(shù)據(jù)描述地形結(jié)構(gòu)[3]。

3）等高線型。等高線通常被存成一個(gè)有序的坐標(biāo)點(diǎn)對(duì)序列，是地形表達(dá)最為常用的形式，能較為直觀科學(xué)地反映地面高程、山體、坡度、坡形、山脈走向等基本的地貌形態(tài)及其變化。利用等高線重新構(gòu)建地表模型時(shí)，只能近似表示真實(shí)的地球表面，在水平切割過程中丟失的詳細(xì)地表信息是不可能從等高線中恢復(fù)的。同時(shí)，單條等高線無法直接反映地貌形態(tài)，必須通過等高線組間接地表示地貌形態(tài)。

4）不規(guī)則三角網(wǎng)（TIN）型。為克服規(guī)則格網(wǎng)的缺點(diǎn)，采用附加地形特征數(shù)據(jù)的方法，按一定的規(guī)則將離散點(diǎn)連接成覆蓋整個(gè)區(qū)域且互不重疊、結(jié)構(gòu)最佳的三角形，從而構(gòu)成完整的DEM，這就是不規(guī)則三角網(wǎng)數(shù)字高程模型[4]。TIN克服了高程矩陣中冗余數(shù)據(jù)的問題，其最主要的優(yōu)點(diǎn)就是具有可變的分辨率，可根據(jù)不同地形，選取合適的采樣點(diǎn)數(shù)。另外，TIN還具有考慮重要表面數(shù)據(jù)點(diǎn)的能力，能利用地貌的特征點(diǎn)線，較好地表示復(fù)雜地形。但其缺點(diǎn)是數(shù)據(jù)結(jié)構(gòu)復(fù)雜，構(gòu)成TIN算法復(fù)雜，計(jì)算時(shí)間長，大區(qū)域的TIN的管理十分困難。

DEM具有表達(dá)樣式多、精度恒定、自動(dòng)化、實(shí)時(shí)性強(qiáng)、尺度綜合性好等特點(diǎn)，在實(shí)際研究與應(yīng)用中，應(yīng)根據(jù)區(qū)域的研究目的和大小，對(duì)數(shù)據(jù)存儲(chǔ)的結(jié)構(gòu)格式進(jìn)行轉(zhuǎn)換，以達(dá)到最優(yōu)化的存儲(chǔ)和快速查詢與編輯。

2 基于DEM的地形要素提取實(shí)驗(yàn)

2.1 實(shí)驗(yàn)區(qū)概況

武夷山地處中國福建省的西北部，江西省東部，位于福建與江西的交界處。地理坐標(biāo)為：北緯27°32'36"～27°55'15"，東經(jīng) 117°24'12"～118°02'50"，總面積999.75 km2，核心面積63 575 hm2，同時(shí)劃定了外圍保護(hù)緩沖區(qū)地帶，面積27 888 hm2。武夷山風(fēng)景名勝區(qū)主要景區(qū)方圓70 km2，平均海拔350 m，屬典型的丹霞地貌，素有“碧水丹山”、“奇秀甲東南”之美譽(yù)。

2.2 DEM應(yīng)用分析

本實(shí)驗(yàn)使用軟件為美國ESRI公司的ArcGIS9.X。首先對(duì)實(shí)驗(yàn)區(qū)元數(shù)據(jù)進(jìn)行瀏覽，得知該數(shù)據(jù)覆蓋范圍，投影為橫軸墨卡托投影，單位為m，其中央經(jīng)線為東經(jīng)117°；再查看圖層屬性信息，在Source選項(xiàng)卡中，可看到該柵格圖層的相關(guān)信息，如該圖層有行列值，柵格大小約為26.17 m。柵格圖層的行列數(shù)和柵格大小是2個(gè)很重要的信息，對(duì)后面的參數(shù)設(shè)定有最直接的影響。

2.2.1 地形指標(biāo)的提取

地形指標(biāo)是描述和評(píng)價(jià)一個(gè)地區(qū)地形情況的有效工具，是進(jìn)行多要素疊加分析中一個(gè)可選條件，也是最短路徑分析中獲取成本圖層的一個(gè)可行方法。地形指標(biāo)一般包括坡度、坡度變率、坡向、坡向變率、起伏度和粗糙度等。

通過ArcMap的“Aspect”工具可先提取出原DEM的坡向信息（見圖1）。而提取坡向變率（SOA）的過程較為復(fù)雜，首先要進(jìn)行反地形的提取，然后對(duì)反地形（見圖2）提取坡向信息，再分別對(duì)正、反地形的坡向進(jìn)行坡度提取（所得結(jié)果為“SOA1”和“SOA2”），最后利用Raster Calculator輸入模型表達(dá)式：SOA=(([SOA1]+ [SOA2]) - Abs([SOA1]- [SOA2]))/ 2，從而得到誤差糾正后的“SOA”。這里需要注意的是，在提取反地形坡向變率時(shí)，需用一個(gè)常數(shù)值減去原DEM，從而得到反地形數(shù)據(jù)。由于求算坡向數(shù)據(jù)與反地形的絕對(duì)高程無關(guān)，僅與相對(duì)高程有關(guān)，因此這里的常數(shù)可任意設(shè)置（一般取原DEM的最大高程值）。

圖1 實(shí)驗(yàn)區(qū)坡向圖層

圖2 實(shí)驗(yàn)區(qū)反地形圖層

2.2.2 等高線的提取

一般來說，在提取等高線時(shí)等高距設(shè)置為大于原DEM柵格大小減30%的區(qū)間為好，這樣可有效避免出現(xiàn)孤立的等高線小尖角、等高線小圓圈等情況。例如，在本實(shí)驗(yàn)中原DEM的柵格大小約為26.17 m，將等高距設(shè)置為大于20 m的數(shù)值，得到的結(jié)果比較符合實(shí)際情況。這也從另一個(gè)角度說明了查看原數(shù)據(jù)信息的重要性（見圖3、圖4）。

圖3 提取后查看100 m等高距

2.2.3 山脊線、山谷線的提取

山脊線和山谷線是地形的重要特征。從地貌學(xué)的角度，它們是坡向變率較大的點(diǎn)的連線；從水文學(xué)的角度，它們則分別是分水線和匯水線[5]。相應(yīng)地，提取這2種地形特征的方法就有空間分析和水文分析2種。

1）空間分析法。先用鄰域統(tǒng)計(jì)分析工具，求算出鄰域平均值圖層；然后用原DEM減去鄰域平均值圖層，得到差值圖層[c]；再通過柵格計(jì)算器，分別求出“sj1=[cha]> 0 & [SOA]> 70，sj2=[cha]> 0 & [SOA]>60”，以及sj3=[cha]> 20 & [SOA]> 60的邏輯結(jié)果。

比較“sj1”、“sj2”和“sj3”3個(gè)圖層，可知設(shè)定不同的閾值，所得到的結(jié)果也是不同的。一般來說，閾值范圍越大，山脊線就越粗越密集，偽山脊線存在的可能性也越大；反之，則山脊線就越細(xì)越稀疏，偽山脊線存在的可能性也就越小。根據(jù)需要和具體情況對(duì)閾值進(jìn)行設(shè)定，從而提取適當(dāng)密度的山脊線和山谷線（見圖5、圖6）。

圖5 提取sj1計(jì)算結(jié)果

圖6 提取sj2計(jì)算結(jié)果

2）水文分析法。先對(duì)原DEM進(jìn)行流向的提取，并利用ArcToolbox中水文分析Sink工具，對(duì)該流向數(shù)據(jù)進(jìn)行洼地的檢查。如存在洼地，則用Fill工具進(jìn)行洼地填充，以消除錯(cuò)誤的流向。對(duì)洼地填充后的結(jié)果進(jìn)行流向提取，并利用流向數(shù)據(jù)進(jìn)行匯流累積量提取，運(yùn)用Raster Calculator提取匯流累積量為0的區(qū)域，并對(duì)其進(jìn)行鄰域平均值平滑。然后，對(duì)平滑后的結(jié)果進(jìn)行重分類。經(jīng)反復(fù)實(shí)驗(yàn)，0.555 4為最佳分界閾值，接近1的部分賦值為1，接近0的部分賦值為0。最后，計(jì)算正地形和前面重分類結(jié)果的乘積，將所得結(jié)果為0的部分重分類為Nodata,1的部分保留為1，從而完成了水文分析法的山脊線提取。山谷線的提取方法與山脊線基本類似，只是用反地形數(shù)據(jù)當(dāng)作原數(shù)據(jù)，重復(fù)山脊線提取步驟，則可得到反地形的山脊線，即原地形的山谷線數(shù)據(jù)（見圖7、圖8）。

圖7 重分類后提取山脊線圖層

圖8 重分類后提取山谷線圖層

3 結(jié) 語

1）對(duì)于相似的要素提取過程或應(yīng)用模型，應(yīng)當(dāng)掌握其共同的算法思想[6,7]，同時(shí)注意區(qū)別其不同之處。在具體應(yīng)用中，需清楚其原理，了解更改每一步的參數(shù)設(shè)置之后會(huì)得到哪些相應(yīng)的變化，對(duì)所得到的結(jié)果應(yīng)結(jié)合已有的數(shù)據(jù)進(jìn)行正誤的判斷和分析。

2）對(duì)于數(shù)據(jù)量較大的工程，在計(jì)算機(jī)設(shè)備有限的條件下，應(yīng)當(dāng)靈活設(shè)置鄰域窗口、輸出柵格尺寸等參數(shù)，在運(yùn)算時(shí)間和結(jié)果精度之間尋找平衡點(diǎn)。

3）對(duì)于水文分析，應(yīng)先對(duì)原DEM的流向數(shù)據(jù)作“Sink”檢驗(yàn)，發(fā)現(xiàn)有洼地存在時(shí)，應(yīng)及時(shí)填充，以消除不合理的河流方向。

4）閾值的設(shè)置是實(shí)驗(yàn)結(jié)果好壞的關(guān)鍵，有些可借鑒其他研究者的經(jīng)驗(yàn)數(shù)值，有些則需實(shí)驗(yàn)比較，從而得出最優(yōu)的閾值。

[1]李志林,朱慶.數(shù)字高程模型[M]. 武漢：武漢大學(xué)出版社,2008

[2]湯國安,劉學(xué)軍,閭國年.數(shù)字高程模型及地學(xué)分析的原理與方法[M].北京:科學(xué)出版社,2005

[3]湯國安,楊昕.ArcGIS地理信息系統(tǒng)空間分析實(shí)驗(yàn)教程[M].北京:科學(xué)出版社,2008

[4]艾廷華,祝國瑞,張根壽.基于Delaunay三角網(wǎng)模型的等高線地形特征提取及谷地樹結(jié)構(gòu)化組織[J].遙感學(xué)報(bào),2004,31(2):43-47

[5]靳海亮,康建榮,高井祥.利用等高線數(shù)據(jù)提取山脊(谷)線算法研究[J].武漢大學(xué)學(xué)報(bào)：信息科學(xué)版,2005(9):810-812

[6]朱慶,趙杰,鐘正.基于規(guī)則格網(wǎng)DEM的地形特征線的提取算法[J].測(cè)繪學(xué)報(bào).2004,33(1):77-82

[7]劉學(xué)軍.基于規(guī)則網(wǎng)數(shù)字高程模型解譯算法誤差分析與評(píng)價(jià)[D].武漢：武漢大學(xué)，2002