基于DEM地形要素提取方法的實驗

孫嘉駿,黃天勇,陳 科

（1.中國水電顧問集團 昆明勘測設計研究院有限公司 測繪地理信息分院,云南 昆明 650041）

DEM是通過有限的地形高程數據實現對地形曲面的數字化模擬或是地形表面形態的數字化表示的量化模型[1]。近年來，隨著空間數據基礎設施的建設和“數字地球”戰略的實施，加快了DEM與地理信息系統、遙感等的一體化進程，為DEM的應用開辟了更廣闊的天地。

1 DEM概述

1.1 DEM的含義

DEM表示三維向量有限序列，是對地球表面地形地貌的一種離散的數字表達，也是測繪學中最常用的地形表示模型，用函數形式描述為：iiii

式中，（Xi,Yi）是平面坐標；Zi是（Xi,Yi）對應的高程。當該序列中各平面向量的平面位置呈規則格網排列時，其平面坐標可省略，此時DEM就簡化為一維向量序列{Zi, i=1,2,3,…,n}[1]。

1.2 DEM的表示模型及特點

隨著對DEM需求的不斷增加，對DEM數據的多樣性、實用性及精度提出了更高的要求。同時由于其數據來源不同，導致數據存儲結構的差別。DEM數據結構表示模型主要包括：

1）不規則離散點型。由不規則離散點表示的DEM是將連續地球表面形態離散成在某一個區域D上的以Xi、Yi、Zi三維坐標形式存儲的高程點Zi((Xi,Yi)ID)的集合[2]。離散點DEM往往是通過測量直接獲取的地球表面的原始特征點數據。這些特征點之間相互獨立，彼此間沒有任何聯系，它是DEM中最簡單的數據組織形式。

2）規則格網型。規則格網DEM是指在水平方向和垂直方向按等間隔方式記錄格網點的坐標，格網點的高程Z表示地形，格網點的平面坐標隱含在行列號中，故適宜用矩陣形式進行存儲，即按行(或列)逐一記錄每一個格網單元的高程值。通常，為了壓縮柵格DEM的冗余數據，可采取用游程編碼或四叉樹編碼方法對數據進行處理；為準確表示地形細部，可采用附加地形特征數據描述地形結構[3]。

3）等高線型。等高線通常被存成一個有序的坐標點對序列，是地形表達最為常用的形式，能較為直觀科學地反映地面高程、山體、坡度、坡形、山脈走向等基本的地貌形態及其變化。利用等高線重新構建地表模型時，只能近似表示真實的地球表面，在水平切割過程中丟失的詳細地表信息是不可能從等高線中恢復的。同時，單條等高線無法直接反映地貌形態，必須通過等高線組間接地表示地貌形態。

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

DEM具有表達樣式多、精度恒定、自動化、實時性強、尺度綜合性好等特點，在實際研究與應用中，應根據區域的研究目的和大小，對數據存儲的結構格式進行轉換，以達到最優化的存儲和快速查詢與編輯。

2 基于DEM的地形要素提取實驗

2.1 實驗區概況

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

2.2 DEM應用分析

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

2.2.1 地形指標的提取

地形指標是描述和評價一個地區地形情況的有效工具，是進行多要素疊加分析中一個可選條件，也是最短路徑分析中獲取成本圖層的一個可行方法。地形指標一般包括坡度、坡度變率、坡向、坡向變率、起伏度和粗糙度等。

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

圖1 實驗區坡向圖層

圖2 實驗區反地形圖層

2.2.2 等高線的提取

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

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

2.2.3 山脊線、山谷線的提取

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

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

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

圖5 提取sj1計算結果

圖6 提取sj2計算結果

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

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

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

3 結 語

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

2）對于數據量較大的工程，在計算機設備有限的條件下，應當靈活設置鄰域窗口、輸出柵格尺寸等參數，在運算時間和結果精度之間尋找平衡點。

3）對于水文分析，應先對原DEM的流向數據作“Sink”檢驗，發現有洼地存在時，應及時填充，以消除不合理的河流方向。

4）閾值的設置是實驗結果好壞的關鍵，有些可借鑒其他研究者的經驗數值，有些則需實驗比較，從而得出最優的閾值。

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

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

[3]湯國安,楊昕.ArcGIS地理信息系統空間分析實驗教程[M].北京:科學出版社,2008

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

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

[6]朱慶,趙杰,鐘正.基于規則格網DEM的地形特征線的提取算法[J].測繪學報.2004,33(1):77-82

[7]劉學軍.基于規則網數字高程模型解譯算法誤差分析與評價[D].武漢：武漢大學，2002