王長海，肖亮亮

廣西交通設(shè)計集團有限公司，廣西南寧 530029

數(shù)字高程模型(digital elevation model, DEM)中包涵海量和多樣的地形特征與構(gòu)造信息，是對地形地貌特征和水文地質(zhì)等信息的空間變化進行定量描述的基礎(chǔ)數(shù)據(jù)[1]，也是自動劃分地形地貌類型的數(shù)據(jù)源. 以DEM數(shù)據(jù)為基礎(chǔ)，以數(shù)字地形分析 (digital terrain analysis, DTA)等相關(guān)的數(shù)字技術(shù)為手段，通過提取分析各類地形地貌特征屬性，自動劃分特征地貌的各種類型，對地貌的自動分類研究有重要價值.

過去人們常常利用目視解譯和人工判讀輔助野外實地勘察實現(xiàn)地形地貌單元的分類， 基于手工目測勾繪出地貌的實體界線. 該方法精度較高，但耗時長且工作量大，無法快速提取大范圍的地貌信息. 此外，還受制圖人員主觀因素影響，人工制作的地貌界線難以保持一致[2]. 在數(shù)字地形技術(shù)快速發(fā)展的前提下，人們已可通過計算機等工具實現(xiàn)地貌單元的自動或半自動劃分. 對象和尺度是地形分類中最重要的兩個因素[3]，然而基于DEM以對象為單位實現(xiàn)對山體對象的提取仍處于研究之中[4].

傳統(tǒng)地貌類型分類方法過分強調(diào)地貌單元的完整性，且無有效的地形地貌實體單元自動提取方法. 本研究針對這個難題，提出了一種基于DEM的山地、平原和嵌套山體界線的自動提取方法，利用坡度等地形因子以及水文分析等相關(guān)方法，能夠?qū)崿F(xiàn)嵌套山體界線和山地、平原的自動化提取.

考慮地形地貌環(huán)境是公路選線設(shè)計的重要原則之一[5]. 將山體提取結(jié)果應(yīng)用到公路選線中，通過將各種地形影響因子轉(zhuǎn)化為成本數(shù)據(jù)集，可以快速、自動地選出最優(yōu)路徑[6]. 相比于僅考慮微觀地貌因子，同時顧及地貌實體單元完整性的選線結(jié)果更加符合實際，這對于減少公路初步選線工作量、優(yōu)化公路選線方法具有重要意義.

1 地形形態(tài)劃分

如何劃分地貌形態(tài)的類型直接關(guān)系到地貌生產(chǎn)利用及地貌形態(tài)的認知科普等問題[7]. 人們通?；诰植康匦蔚某ＷR來進行地形地貌分類，通過大量從定性到定量的研究形成了多個地形分類體系和標(biāo)準(zhǔn)，主要包括基于地貌學(xué)、水文學(xué)、自然地理學(xué)和土壤學(xué)等分類體系.

1.1 地形形態(tài)劃分體系

地形地貌的分類始于人們對局部地形的主觀認識，主要包括基于水文學(xué)、自然地理學(xué)、土壤學(xué)及地貌學(xué)等分類體系.

基于水文學(xué)的分類體系主要根據(jù)地形所受的水文沖刷侵蝕的不同程度，基于縱坡面的坡長、坡寬及坡度，將地形坡面分為山肩、山頂、背坡、趾坡、麓坡和沖擊地等.

根據(jù)自然地理學(xué)，地表形態(tài)差異必然引起各種自然地理現(xiàn)象及過程的變化[8]. 宏觀上，地形可劃分為平原、丘陵、高山、低山和高原等類型； 微觀上，斜坡可劃分為緩坡、中坡和陡坡等. 這種劃分體系完全是以坡度為劃分指標(biāo)，常見的有以15°、25°和45°為臨界值對坡體進行連續(xù)劃分. 因為坡度是地形表面的單點屬性之一，所以這種分類體系對于利用柵格DEM實現(xiàn)地貌單元的劃分十分有利，本研究對平原、山地的劃分也利用到了這一劃分體系.

基于地貌形態(tài)學(xué)，將地形形態(tài)特征作為研究對象，依據(jù)描述要素的空間組合與分析比對來進行特征地貌的獲取[9]. 基于這種方案，任何復(fù)雜的地形地貌均可以劃分成3種類型的要素，分別為平地、洼地和坡地.

1.2 地形形態(tài)分類指標(biāo)

利用DTA技術(shù)基于DEM來實現(xiàn)地形地貌的自動分類必須要給計算機確定的解析規(guī)則和規(guī)范，也就是必須有確定的分類指標(biāo)，即用于描述地貌的地形因子. 地形因子是地形信息的載體和最重要的表現(xiàn)形式，也是地形地貌在某一方面的特征的數(shù)字描述. 因此，各種地形因子的準(zhǔn)確提取對地貌學(xué)研究具有重要價值. 地形因子通常可以分為宏觀地形因子和微觀地形因子，如圖1所示.

圖1 地形因子分類框架Fig.1 Classification structure of terrain factor

1.3 地形形態(tài)分類方法

地貌類型劃分時，將人為經(jīng)驗與地貌分類的量化指標(biāo)進行結(jié)合，從而判斷出地貌的實體單元類型，這是人工判讀分類中的常用方法. 此方法在分類時人為的主觀因素仍然占有很大部分，無法進行嚴密而又精確的地形分類. 黃杏元[10]提出了地形分類決策表方法，首先基于區(qū)域地形地貌的特點形成地形分類決策表，獲得分類時所需的地形要素，然后基于DEM進行地形要素提取，再進行地形類型的自動提取，即可獲得區(qū)域的地貌類型，如圖2所示. 該方法為基于DEM自動提取地貌類型提供了最基礎(chǔ)的方法，并且分類過程簡單、快速. 高玄彧[11]提出一種以相對高度為主，絕對高度為輔的地形分類方法. 這種分類方法避免了其他一些分類方法中與人們的觀念相悖的低山不低、高山不高的現(xiàn)象，符合人們對于高山、中山和低山等的觀念.

圖2 自動提取地形信息[10]Fig.2 Automatic extraction of terrain infomation[10]

本研究認為，基于柵格DEM利用DTA技術(shù)對地貌形態(tài)進行劃分的核心問題，在于如何準(zhǔn)確并且高效地從DEM中提取地形因子，如何利用地形因子準(zhǔn)確表達各種地形地貌. 不同的地形形態(tài)分類方法有各自的側(cè)重點，可以根據(jù)具體的應(yīng)用目的選擇不同的分類方法. 但是不管選用什么分類方法，最為重要的還是準(zhǔn)確計算地表形態(tài)要素和擬定有效的地形分類決策方案.

2 地形因子的計算

2.1 坡度和坡向的計算

坡度(slope)和坡向是描述地形復(fù)雜度的重要指標(biāo)，并且相互聯(lián)系. 坡度是地表高度變化率的度量，直接反映斜坡的傾斜程度. 坡向則是斜坡方向的度量，用于反映斜坡所面對的方向[12]. 坡度是地表在某個點的傾斜程度，是經(jīng)過該點的地表微分單元的法矢量n與z軸的夾角，即

(1)

在進行坡度的實際計算時，一般采用簡化的差分公式，即

(2)

其中，fx是x方向的高程變化率；fy是y方向的高程變化率.

地表上的任意一點均有坡度，其單位一般為“度”或“坡度百分數(shù)”. 由坡度計算示意圖(圖3)可知，坡度百分數(shù)為(h/d)×100%， 度為arctan(h/d). 斜坡所面對的方向是坡向，通過某一點法線正方向的平面投影和正北方向之間的夾角則為該點的坡向[13]. 坡向α的數(shù)學(xué)表示為

(3)

圖3 坡度計算示意圖Fig.3 Slope calculation schematic

差分坡度和坡向的計算原理見圖4，設(shè)中心格網(wǎng)點(i,j)的坐標(biāo)為(xi,yj)， 局部地形曲面設(shè)為z=f(x,y)， 格網(wǎng)間距為g.

圖4 差分坡度和坡向計算原理Fig.4 The calculation principle of differential slope and aspect

則有二階差分公式為

(4)

2.2 水文因子的計算

水文因子指在流域分析或者水文分析中所涉及到的一些與地形有關(guān)的因子，如水流長度、水流方向、河流網(wǎng)絡(luò)和匯流累積量等.

洼地是DEM局部地形單元， 被較高的高程所包圍. 在DEM中的洼地是由數(shù)據(jù)采集和內(nèi)插時的誤差造成的，稱為“偽洼地”. 在進行數(shù)字地形的分析時，地形特征線提取的正確性將受偽洼地的影響. 因此，在進行水文因子分析之前需要對DEM進行無洼地化處理，盡可能地減少“偽洼地”的影響，確保所提取的地形特征線能夠保持良好的連續(xù)性，且無離散的點及局部的短線. 一般通過將洼地邊緣的最小高程值賦以洼地區(qū)域中的每一個格網(wǎng)點，從而達到填平洼地的目的.

基于格網(wǎng)DEM的水流分析，其水流方向就是水流離開格網(wǎng)時的流向. 在流域分析中，八方向法是當(dāng)前被廣泛使用的確定流向的方法，該方法將像元的流向指向其周圍8個像元中距離權(quán)重坡度最大的像元(圖5). 該法的局限性在于不允許水流分散到多個像元[14].

圖5 水流方向計算Fig.5 The calculation of flow direction

通過計算中心像元與其8個鄰域像元的距離權(quán)重坡度，可以確定圖5(a)的中心像元流向. 對于任意4個緊鄰的像元，把其中心像元與鄰接點的高程差除以1，可計算得到其坡度. 對于任意4個角落的鄰接像元，其坡度計算則是用高程除以1.414，如圖5(b). 結(jié)果顯示，最陡的坡度是從中心像元到其右側(cè)像元，因此為流向. 八方向法在回流區(qū)和邊界明確的峽谷地區(qū)能夠獲得較好的結(jié)果[15].

圖6 匯流累積量計算Fig.6 The calculation of confluence accumulation

區(qū)域地表的每點流水累積量用匯流累積量數(shù)矩陣來表示. 在進行地表徑流的模擬時， 通過水流方向可以計算得到匯流累積量. 匯流累積量的計算過程如圖6. 圖6(a)為已填洼的DEM，圖6(b)為水流方向柵格，圖6(c)為匯流累積量柵格. 圖6(c)中2個陰影柵格的流量累積值同為2. 其中，上面的柵格像元接受來自它左邊和左下方像元的水流，下面的像元接受來自其左下方像元的水流，且其本身已有流量累計值1. 利用水流方向柵格，匯流累積量柵格列出流經(jīng)它的像元數(shù)，每個格網(wǎng)記錄了有多少個上游的像元將水排給它(不包括被計算的柵格像元本身).

通過匯流累積量柵格可以用兩種方式來解釋：① 匯流累積量高的像元，一般對應(yīng)河道，而山脊線的像元累積量通常為0；② 若乘以每個像元的大小，則其像元值就是排水面積.

3 實驗研究

3.1 實驗數(shù)據(jù)介紹

本研究利用ArcGIS實現(xiàn)對平原、山地高效準(zhǔn)確的提取. 由于四川盆地的地形特征整體地形四周高、中間低，呈盆地地形，西部為高原山地，東部為盆地. 盆地內(nèi)的地形主要是丘陵和平原，平原和山地的地形特征較為顯著. 因此，本實驗采用的數(shù)據(jù)源是四川某地區(qū)的90 m和30 m分辨率的柵格DEM.

3.2 平原/山地自動劃分

為了實驗需要，裁減了部分比較典型的山地和平原的地形數(shù)據(jù)，原始數(shù)據(jù)是img格式的柵格DEM，具體實驗數(shù)據(jù)如圖7所示.

圖7 實驗樣區(qū)數(shù)據(jù)Fig.7 (Color online) Experimental data

根據(jù)《中國1∶100萬地貌制圖規(guī)范》[16]， 丘陵和山地指的是圖斑最高點與圖斑邊緣最高點的高差大于30 m坡度大于7°(一般大于10°)的形態(tài). 根據(jù)平原的構(gòu)成形態(tài)，平原可抽象為多個坡度在一定閾值內(nèi)、相互連接且達到一定面積的坡面所共同組成. 本研究同時對坡度、面積兩個指標(biāo)進行綜合考慮，通過設(shè)置恰當(dāng)?shù)拈撝祦磉M行山地和平原的區(qū)分. 為此，本研究還設(shè)計了相應(yīng)的計算流程來進行平原/山地的自動劃分.

通過多次試驗與反復(fù)對比后，本研究最后確定的坡度閾值為8°，即將原始坡度柵格重分類為0°～8°和8°～90°兩類，此時可以較好地判別山體與平原的界線. 然后通過統(tǒng)計各個區(qū)域的面積，并確定面積閾值. 根據(jù)制圖要求，在1∶100萬比例尺下成圖最小的圖斑所對應(yīng)的實際面積約為4 km2，但考慮本研究地區(qū)地形地貌的結(jié)構(gòu)特點，并經(jīng)反復(fù)實驗及對比目視結(jié)果，最后確定本研究區(qū)域的面積閾值為8 km2. 由于自動提取的結(jié)果與實際地形有一定差異，在自動提取地貌單元之后，還要結(jié)合專家知識來進行制圖綜合，最終劃分結(jié)果如圖8所示.

圖8 山體、平原劃分結(jié)果Fig.8 (Color online) The classification results of mountain and plain

3.3 嵌套山體界線提取

一般情況下，平原和山地的劃分并不能提取山體之間的嵌套關(guān)系，因此，需要對嵌套山體界線做進一步的劃分. 根據(jù)制圖規(guī)范，若兩個相連的山體單元兩側(cè)圖斑與埡口的最大海拔高差大于500 m，原則上應(yīng)將該山體單位劃分為兩個實體. 根據(jù)些規(guī)則對山體進行劃分仍需要相應(yīng)的先驗知識，但它為我們提供了劃分嵌套山體界線的思路.

分界點位置的計算是山體界線提取的關(guān)鍵，得到分界點的位置后，各分界點之間的連線即為山體的分界線. 如圖9(a)所示為相鄰的山體剖面圖. B和D是埡口點，A、C和E是各山峰頂點. A到B的垂直距離為a， B到C的垂直距離為b， C到D的垂直距離為c， D到E的垂直距離為d. 假設(shè)a>500 m、b>500 m、c<500 m和d<500 m. 根據(jù)前面的定義，如果埡口兩側(cè)的最高海拔高差大于500 m，那么此處應(yīng)該為山體的分界點. 但是B點和D點的情況則不同，兩側(cè)的山峰頂點到最低點埡口的距離一側(cè)大于500 m，一側(cè)小于500 m，根據(jù)上述規(guī)則此處不是山體分界點. 所以，根據(jù)a、b、c和d的長度及其相互關(guān)系，可以判定某一埡口點是否為山體的分界點，進而可以將山體分界點的埡口點的求取，轉(zhuǎn)化為求取a、b、c和d的長度[17].

圖9 山體分界點判別方法示意圖[17]Fig.9 Schematic of the cut-off point discrimination method[17]

然而，實際的山地情況比上述情況要復(fù)雜的多，在三維空間范圍內(nèi)，山體在各個方向都有相鄰的連接的山地，因而求解山地分界點的過程變得十分復(fù)雜. 若把山體倒轉(zhuǎn)，則需要將進行判別的埡口點轉(zhuǎn)化成倒轉(zhuǎn)后的山頂點，但其判別的規(guī)則不變. 如圖9(b)所示，該地形為圖9(a)中倒轉(zhuǎn)的地形. 從這個角度看，可以將問題轉(zhuǎn)換為對直觀水流現(xiàn)象的分析. 如圖9(b)所示，山體倒轉(zhuǎn)后將會成為盆地. 由于在a、b、c、d中d最小，如果不斷有水從上方流入該區(qū)域，隨著水量增多，當(dāng)水量超過盆地的容量后，即淹沒深度超過d以后，E點出盆地中的水就會流向C點;同理，由于a、b均大于500 m，故當(dāng)淹沒深度達到500 m時，A、C兩處盆地中的水流相互獨立且互不連通，則B為山體的一個分界點.

因此，根據(jù)以上描述，通過將地形倒轉(zhuǎn)，借助水文淹沒分析，可以將判別規(guī)則從二維情況推廣到三維情況，再結(jié)合匯流累計量等水文因子，即可實現(xiàn)嵌套山體界線的提取. 具體實驗區(qū)域的數(shù)據(jù)展示如圖10所示.

圖10 實驗區(qū)域數(shù)據(jù)Fig.10 Experimental data

圖11 不同閾值所對應(yīng)的水流路徑Fig.11 (Color online) Flow path corresponding to different thresholds

基于倒轉(zhuǎn)地形的水文淹沒分析，將最大海拔高差低于500 m的盆地都填平，依據(jù)山體判別規(guī)則，可以消除掉不滿足的埡口點. 再由匯流累積量得到水流路徑，然后通過設(shè)置不同的匯流累積量閾值得到不同級別的水流路徑. 不同級別的水流路徑對應(yīng)不同的閾值，不同級別的水流路徑作為不同等級的嵌套山體的劃分界線. 不同閾值提取的水流路徑如圖11所示. 實驗結(jié)果表明，通過設(shè)置合理的匯流累計量閾值，可以得到較好的嵌套山體界線的劃分結(jié)果，分界線較好地反映了山體的結(jié)構(gòu)特征.

由圖12可見，采用本研究方法提取的山體界線與依據(jù)專家知識手工勾繪的山體界線吻合較好. 當(dāng)閾值設(shè)置為5 000時，根據(jù)水流路徑和圖廓邊界可以得到最高級別的山體界線，可以將實驗區(qū)域的山體劃分為5個部分. 當(dāng)閾值設(shè)置為2 000時，又可以將該區(qū)域的山體再次劃分為4個子區(qū)域山體，依此類推，根據(jù)其他閾值得到的水流路徑可以得到其他較低級別的山體界線，由此便實現(xiàn)了嵌套山體界線的提取，如圖12所示.

圖12 嵌套山體界線Fig.12 (Color online) Nested mountain boundaries

在提取嵌套山體界線時，依據(jù)上述方法，較好地實現(xiàn)了第1等級和第2等級的山體界線的提取，但對于更低等級的嵌套山體界線，由于基于水流路徑的提取存在某些界線無法閉合的問題，因此仍需綜合利用專家知識、地貌特點和參考特征線等信息，進一步調(diào)整和修改自動提取的界線.

3.4 公路選線

利用ArcGIS的空間分析功能，基于DEM數(shù)據(jù)，綜合考慮地形、坡度和河網(wǎng)分布等各類地形因子信息，再分別按照對地形的不同影響，將其重新分為10個等級. 同時，將實驗區(qū)域整體的山體與平原的分布情況考慮在內(nèi)，利用平原與山體劃分實驗結(jié)果，連同將上述微觀地形因子數(shù)據(jù)集作為公路選線的綜合成本數(shù)據(jù)集. 根據(jù)文獻[6]，坡度、起伏度及河網(wǎng)分布對公路選線影響的重要性依次為“坡度>起伏度>河網(wǎng)分布”. 由于宏觀的山體、平原分布情況，是從宏觀角度考慮地形起伏對公路選線的影響，相比微觀因子，其在實際公路選線中作為預(yù)先初步判斷條件，因此將其實驗重要性排在河網(wǎng)分布之后. 在此基礎(chǔ)上，選定起點和終點，利用ArcGIS計算得到距離成本數(shù)據(jù)集與方向數(shù)據(jù)集，如圖13所示. 最后，再執(zhí)行最短路徑功能得到公路選線的最優(yōu)路徑. 實驗結(jié)果如圖14所示.

圖13 距離成本數(shù)據(jù)集Fig.13 (Color online) Distance costs dataset

圖14 公路選線結(jié)果Fig.14 (Color online) Highway route selection results

根據(jù)實驗結(jié)果，利用地形分析的手段進行公路選線，在計算成本數(shù)據(jù)集時，顧及宏觀山體分布因子與否得到的公路選線結(jié)果出現(xiàn)了一定的偏差. 如圖14所示，從選線結(jié)果路線長度來看，估計宏觀山體因子得到的選線結(jié)果，路線長度比為顧及宏觀山體因子得到的路線更短，并且路線更為筆直，彎曲更少. 從路線經(jīng)過的地形看，顧及宏觀山體因子得到的結(jié)果所經(jīng)過的地形相對更趨一致，路線連續(xù)的起伏相對較少. 總體來看，相比未考慮宏觀山體因子的結(jié)果，估計宏觀山體因子的公路選線結(jié)果更為合理.

結(jié) 論

本研究基于DEM數(shù)據(jù)，通過數(shù)字地形分析方法，利用面積和坡度兩個指標(biāo)，計算出大于一定閾值的坡度，且相互連通的柵格的面積，實現(xiàn)了山地和平原的自動劃分. 研究利用水文淹沒分析法消除了不符合山體判別規(guī)則的山體埡口，且利用水文分析的相關(guān)方法，通過設(shè)置合適的匯流累計量閾值，較好地實現(xiàn)了山體界線的自動劃分，實驗結(jié)果與依據(jù)專家知識手工勾繪的結(jié)果相比，取得了良好的效果. 基于該方法提取的山體界線，結(jié)合專家知識加以修正及制圖綜合，能夠加快地貌制圖的速度，提高地貌基本形態(tài)類型的分類精度，為大面積地貌信息的高效提取與綜合利用奠定基礎(chǔ). 最后，本研究基于山地和平原的劃分結(jié)果，分別在顧及宏觀山體因子與否的情況下進行了公路選線實驗. 從實驗結(jié)果來看，與傳統(tǒng)只考慮微觀地形因子的公路選線結(jié)果相比，顧及宏觀山體因子的公路選線結(jié)果更為合理，可以為初步的公路選線提供更加可靠合理的結(jié)果，大大提高公路選線的工作效率.

本研究所采用的方法在嵌套山體提取時，坡度和匯流累積量閾值的確定以及山體界線的閉合需要一定的人工干預(yù)，存在自動化程度不足的問題. 在以后的研究中擬增加曲率、鞍部點和山脊線等指標(biāo)，擬通過建立地形對象的層次關(guān)系，真正實現(xiàn)以對象為單位的地貌實體界線的自動化提取. 在公路選線中，擬增加地質(zhì)條件、土地利用等因素，同時研究不同因素對公路選線成本的權(quán)重，得到更加合理的公路選線結(jié)果.