基于DEM提取地性線研究

栗 業(yè)

（河北省地質(zhì)礦產(chǎn)勘查開發(fā)局第二地質(zhì)大隊(duì)，河北 唐山 063000）

1 研究背景概述

1.1 研究背景

數(shù)字高程模型（Digital Elevation Model， DEM）是依據(jù)高程數(shù)據(jù)對地形曲面的數(shù)字化模擬，是地表形態(tài)數(shù)字化的量化模型。DEM是地理空間地形分析的核心數(shù)據(jù)，是新一代的地形圖，同時(shí)也是GIS數(shù)據(jù)庫中重要的組成部分。DEM與傳統(tǒng)地形圖相比有著巨大優(yōu)勢。首先，它可以在計(jì)算機(jī)中直接應(yīng)用、計(jì)算，為各種輔助設(shè)計(jì)系統(tǒng)利用；其次，DEM是以數(shù)字形式存儲，便于數(shù)據(jù)的錄入、修改、更新、復(fù)制和管理，同時(shí)，可以方便地轉(zhuǎn)換成其他形式的材料文件和產(chǎn)品；最后，DEM中含有傳統(tǒng)地形圖中所無法容納和表達(dá)的垂直分布地物信息，即高程信息，所以DEM數(shù)據(jù)可以滿足眾多行業(yè)的發(fā)展需求。

1.2 研究意義

描述地貌形態(tài)的骨架線、控制線稱為地性線，包括山脊線、山谷線、鞍部、凹地等線性物體。山谷線和山脊線構(gòu)成了地勢起伏變化的特征線，對地形地貌研究有重要意義。DEM數(shù)字高程模型中表達(dá)地形地貌的特征線，即地性線，其中包含了大量的地性線信息?；贒EM提取地性線的研究在區(qū)域水文分析、土木工程建設(shè)、地質(zhì)勘探選定資源靶區(qū)、遙感影像自動(dòng)配準(zhǔn)等諸多領(lǐng)域具有重要的實(shí)踐應(yīng)用價(jià)值和科學(xué)意義。

2 國內(nèi)外研究現(xiàn)狀

多年來，國內(nèi)外專家學(xué)者在數(shù)字高程模型中提取隱含著的山脊線和山谷線信息，及其如何將其快速自動(dòng)提取，并應(yīng)用于水文分析、地形分析、3D視域等多個(gè)研究領(lǐng)域方面，取得了一定的研究成果。

1984年，Yeoli提出了一個(gè)兩步算法。首先，使用斷面識別方法確定特征明顯的點(diǎn)，存儲于XMIN、YMIN、HMIN這3個(gè)矩陣中；然后從上述3個(gè)矩陣中去除已有山谷線上的點(diǎn)，用剩余的點(diǎn)中最大的點(diǎn)開始一條新山谷線的搜索。1991年，John Fairfield等在地表流水模型的基礎(chǔ)上，推出了隨機(jī)水流模型。該模型基本解決了DEM網(wǎng)格問題，但數(shù)據(jù)量成幾十倍的增加，計(jì)算耗時(shí)多，且只適用于區(qū)域較小且邊界較規(guī)則的范圍；1995年，Amnon Meisels等提出了骨架法。該算法是提取等值線上曲率值較大的點(diǎn)，然后把它們連接起來作為特征線。這種算法所需的參數(shù)少，缺點(diǎn)是不能在較大凹陷地區(qū)使用。

近年來，國內(nèi)學(xué)者也在此領(lǐng)域取得了重要突破。張維軍應(yīng)用道格拉斯-帕克提取山脊和山谷的特征點(diǎn)，并結(jié)合提取區(qū)域地形大致劃分和匯水線的幾何分析方法， 用概略的分水線和匯水線對特征點(diǎn)進(jìn)行識別、歸類、順序提取各條山脊線和山谷線。姚慧敏等闡述了利用數(shù)學(xué)形態(tài)學(xué)的流域分割算法提取地性線的原理， 采用矢量柵格混合算法提取地性線的具體步驟， 最后實(shí)驗(yàn)結(jié)果表明， 其提出的算法是可行的， 提取的地性線有較好的連續(xù)性。黃培之研究了特定實(shí)驗(yàn)區(qū)域的山脊線和山谷線的幾何特性或物理特性的單一方面設(shè)計(jì)的提取山脊線和山谷線的算法，隨后又提出了一種基于地形表面流水分析與等高線幾何分析結(jié)合的提取地性線的方法。該方法把等高線幾何分析的方法與地形表面流水模擬分析的方法有機(jī)地結(jié)合起來， 能夠克服各自所具有的弊端。

3 DEM與地性線的數(shù)學(xué)表達(dá)

3.1 DEM的數(shù)學(xué)表達(dá)及數(shù)據(jù)結(jié)構(gòu)

DEM數(shù)字高程模型是地表形態(tài)的數(shù)字化表達(dá)。從狹義角度說DEM是區(qū)域地表正高的數(shù)字化表達(dá)模型；廣義上來講，它是時(shí)空中的地理對象表面正高高度的數(shù)字化表達(dá)模型。

數(shù)學(xué)意義上的DEM是表示區(qū)域A上的三維向量有限序列，函數(shù)表達(dá)為：

式中，（Xi，Yi）是平面坐標(biāo)，Zi是（Xi，Yi）對應(yīng)的高程。

目前，DEM主要有4種數(shù)據(jù)結(jié)構(gòu)：①不規(guī)則離散點(diǎn)數(shù)據(jù)結(jié)構(gòu)。不規(guī)則離散點(diǎn)是最簡單的組織形式，它由不規(guī)則的離散點(diǎn)表示地表形態(tài)。這些點(diǎn)是通過測量獲得的獨(dú)立原始特征點(diǎn)數(shù)據(jù)。②規(guī)則格網(wǎng)（Grid）數(shù)據(jù)結(jié)構(gòu)。在垂直、水平兩個(gè)方向上以相等間隔逐行逐列記錄每個(gè)網(wǎng)格的高程值稱為規(guī)則格網(wǎng)DEM。為方便行列號和平面坐標(biāo)的轉(zhuǎn)換，規(guī)則網(wǎng)格中還標(biāo)識了西南角坐標(biāo)值、網(wǎng)格間距等。為壓縮存儲空間，常將規(guī)則網(wǎng)格每個(gè)單元值減去區(qū)域內(nèi)最小高程值，再乘以一個(gè)常數(shù)是規(guī)則網(wǎng)格值全變?yōu)檎麛?shù)。規(guī)則格網(wǎng)DEM通常包含記錄西南角起點(diǎn)坐標(biāo)、坐標(biāo)類型、行列數(shù)的頭文件和記錄行列高程矩陣的數(shù)據(jù)體。規(guī)則格網(wǎng)的編碼結(jié)構(gòu)有行程編碼、塊狀編碼、四叉樹編碼。③不規(guī)則三角網(wǎng)（TIN）數(shù)據(jù)結(jié)構(gòu)。不規(guī)則三角網(wǎng)就是一系列互不重疊互不相交的三角形連接在一起表示地表形態(tài)。不規(guī)則三角網(wǎng)一般采用鏈表存儲，由結(jié)點(diǎn)和三角形頂點(diǎn)兩組記錄組成。TIN是典型的矢量拓?fù)浣Y(jié)構(gòu)。TIN克服了規(guī)則格網(wǎng)的缺點(diǎn)，數(shù)據(jù)量精簡，具有可變的分辨率，能表達(dá)復(fù)雜的地形。但該數(shù)據(jù)結(jié)構(gòu)、算法都較復(fù)雜，管理困難，計(jì)算耗時(shí)長，適用于高精度小范圍的地形建模。④等高線數(shù)據(jù)結(jié)構(gòu)。等高線是最常用的地形表達(dá)形式，在地形圖上高程相等的點(diǎn)所連成的閉合曲線即為等高線。它直觀地反映了高程、坡度、坡向、山脈走向等基本地貌。等高線以坐標(biāo)點(diǎn)對序列形式存儲，在利用等高線重建地表模型時(shí)容易丟失細(xì)節(jié)信息，而且單條等高線無法反映地貌。

3.2 地性線的數(shù)學(xué)表達(dá)

3.2.1 山脊線

在實(shí)際地形地貌中條帶狀隆起的頂部線即山脊線。曲面上二階方向?qū)?shù)取得極小值，且一階導(dǎo)數(shù)為0的點(diǎn)是山脊點(diǎn)。山脊線則是相鄰山脊點(diǎn)的連線。

假設(shè)Ta和Tb正負(fù)型反斜坡，Lab=Ta∩Tb是兩斜坡的交線。在交線上存在一點(diǎn)P，其高程為Hp，則對任意iLab，存在CMi=Bi∩Ta∩Tb，使 得Hp-Hi≥0，得 出T=Ta∪Tb為山脊，Lab為山脊線或稱分水線，Bi為過點(diǎn)i的平面，CMi為山脊線等高線。如果地形曲面T的表達(dá)式為Z = A00+A10X+ A01Y+A11XY+A20X2+A02Y2+……（A20<0；A02<0且A11≠0） 時(shí)，判斷點(diǎn)P是否為山脊上一點(diǎn)，需滿足任意下面條件：坡度α=0；最大曲率Cmax>0；且最小曲率Cmin=0或坡度α>0；且斷面曲率Cs>0。

3.2.2 山谷線

山脊線和山谷線是兩個(gè)完全對應(yīng)相反的概念。曲面上二階方向?qū)?shù)取得極大值，且一階導(dǎo)數(shù)為0的點(diǎn)是山谷點(diǎn)。山谷線則是相鄰山脊點(diǎn)的連線。山谷線收集兩側(cè)的水，屬于匯水區(qū)域，易形成河流。山谷線是純粹的局部地貌概念，因?yàn)樗诙它c(diǎn)處的一小段范圍內(nèi)可能不會形成水系。由于它的局部性質(zhì)，它比河流更適合，更便于計(jì)算機(jī)提取。所有我們近似地把它當(dāng)成河流的一部分。

假設(shè)Ta和Tb正負(fù)型反斜坡，Lab=Ta∩Tb是兩斜坡的交線。在交線上存在一點(diǎn)P，其高程為Hp，則對任意iLab，存在CMi=Bi∩Ta∩Tb，使 得Hp-Hi≤0，得 出T=TaTb為 谷底，Lab為山谷線或稱匯水線，Bi為過點(diǎn)i的平面，CMi為山谷線等高線。如果地形曲面T的表達(dá)式為Z=A00+A10X+ A01Y+A11XY+A20X2+A02Y2+……（A20>0； A02>0且A11≠0） 時(shí)，判斷點(diǎn)P是否為山谷上一點(diǎn)，需滿足任意下面條件：坡度α=0；最大曲率Cmax=0；且最小曲率Cmin<0或坡度α>0；且斷面曲率Cs<0。

4 基于DEM提取地性線的數(shù)學(xué)算法

4.1 基于圖像處理的算法

DEM實(shí)質(zhì)是一種柵格數(shù)據(jù)形式，網(wǎng)格大小即為像素大小，高程值為灰度值。DEM數(shù)據(jù)可以按數(shù)字圖像處理的方法來設(shè)計(jì)算子。在DEM中高程值越大，灰度越大，山脊線上的點(diǎn)比其他地方的點(diǎn)都要亮。按照各種濾波算子就能進(jìn)行邊緣提取，將地形特征線提取出來。Chakrevavanieh提出用高斯算子從DEM上自動(dòng)提取山脊線和山谷線上的點(diǎn)，并將其用于規(guī)則格網(wǎng)的數(shù)據(jù)壓縮。唐心紅等在提取地性線前利用小波分析進(jìn)行多尺度表示，在多尺度DEM中提取。余生晨利用傅里葉變換在頻域中計(jì)算山脊線和山谷線，但該方法計(jì)算量大?；趫D像處理的技術(shù)還有以下兩個(gè)方面的難點(diǎn)：排除DEM中的噪聲影響；所提取特征點(diǎn)連接成線的算法問題。

4.2 基于地表幾何形態(tài)的算法

斷面極值法是地表幾何形態(tài)分析的基礎(chǔ)算法。該算法的基本思路是：在DEM垂直和水平橫截面上的極值，極大值就是分水點(diǎn)，極小值是匯水點(diǎn)；然后根據(jù)相應(yīng)的準(zhǔn)則劃分極值點(diǎn)歸屬形成地性線。地表幾何形態(tài)分析算法由于只判斷兩個(gè)斷面，提取過程中會出現(xiàn)遺漏現(xiàn)象。而且閾值過大會丟失特征點(diǎn)，形成的地性線也就斷裂較多，長度較短，閾值過小，特征點(diǎn)會出現(xiàn)過度提取現(xiàn)象，地性線冗余。

Fowler改進(jìn)該算法，提出用4柵格掃面標(biāo)記山脊點(diǎn)山谷點(diǎn)，這4柵格單元分別為（i， j）、（i， j+1）、（i+1， j）、（i+1， j+1），相比于8鄰域算法效率高而且不會遺漏山脊點(diǎn)和山谷點(diǎn)，缺陷是如果DEM采樣點(diǎn)密集會提取多余的特征點(diǎn)。Jenson認(rèn)為任意一柵格與相鄰8單元格形成4對位置，至少一對高程都比該點(diǎn)高程值大，則該點(diǎn)可能是山谷點(diǎn)，如果至少有一對高程點(diǎn)比該點(diǎn)小，該點(diǎn)就可能為山脊點(diǎn)。這種改進(jìn)算法易于編程實(shí)現(xiàn)，但產(chǎn)生地性線的斷裂較多。針對特征點(diǎn)丟失現(xiàn)象，黃培之提出再增加一組對角線上斷面的方法來克服這一缺陷，研究實(shí)驗(yàn)證明該法有效，但同時(shí)增加了特征點(diǎn)噪聲。

4.3 基于地形表面流水物理模型分析算法

根據(jù)水往低處流，最終順谷而流，匯集成水系網(wǎng)的自然規(guī)律，Chakrevavanieh率先提出模擬地表流水物理模型。D8法是該模型中的最常用算法：假設(shè)只有一個(gè)單一的網(wǎng)格的水流進(jìn)八個(gè)相鄰的網(wǎng)格單元。水流方向以最陡坡度法來確定，即在3×3的DEM柵格上，計(jì)算中心柵格與各相鄰柵格間的距離權(quán)落差，取距離權(quán)落差最大的柵格為中心柵格的流出柵格。以數(shù)值表示每個(gè)單元的流向，變化范圍是1~255。然后匯流柵格圖由流向柵格圖生成。匯流柵格上每個(gè)單元的值代表上游匯流區(qū)匯入該單元的流入路徑數(shù)較大者。從流域匯流柵格圖中可以輕易地提取流域各種參數(shù)特征值。在實(shí)際應(yīng)用時(shí)需要解決以下問題：①水道起始點(diǎn)位置的確定；②偽負(fù)地形的識別及填充；③凹陷與平坦處水流方向的設(shè)定。而且基于地表水流模擬方法存在兩個(gè)缺陷：①洼地的填充和平坦地區(qū)水流方向的設(shè)定十分麻煩；②該方法所計(jì)算的匯水量與高程有關(guān)，計(jì)算的結(jié)果為高程值大的區(qū)域匯水量小，高程值小的區(qū)域匯水量大，由于水流遵從高處連續(xù)性的向下流動(dòng)的自然特征，這就使得在一個(gè)河道水系或山谷線上的相對位置較高的某一點(diǎn)因其匯水量較小而被遺漏，而再非山谷處處于低位置的某點(diǎn)因其匯水量較大而被誤判為山谷線上的點(diǎn)。在后來的研究中，許多學(xué)者提出利用平滑DEM，墊高填平凹陷區(qū)，標(biāo)識洼地集水區(qū)的方法解決上述問題，并取得了良好效果。

4.4 基于地形表面幾何形態(tài)分析和地表流水物理模型分析的算法

由于地表幾何形態(tài)分析和地表流水分析模型的算法有一些不足之處，所以專家們提出了將兩種模式結(jié)合提取地形線。這種算法的基本思路是首先按照地表流水模型在稀疏的DEM中概略提取地性線；然后對周圍地域進(jìn)行幾何地理形態(tài)分析，在精確定位地性線。該方法的關(guān)鍵是求出概略地性線與DEM網(wǎng)格的交點(diǎn)，在該點(diǎn)周圍的小范圍內(nèi)進(jìn)行幾何分析，找出正交方向地形斷面上的高程極值點(diǎn)，正是這一點(diǎn)為精確點(diǎn)。這種綜合算法在特征點(diǎn)提取閾值的選取、DEM網(wǎng)格大小的確定方面還有待進(jìn)一步研究。

4.5 平面曲率與坡形組合法

求取DEM平面曲率和地表的正負(fù)地形，正地形上曲率的大值是山脊，負(fù)地形上曲率的大值是山谷。這種方法中平面曲率的大小可以用來調(diào)教提取寬度，方法簡便效果好，但提取的圖形是柵格形式，在轉(zhuǎn)為矢量線性方面還需進(jìn)一步研究改進(jìn)。

5 結(jié)語

綜上所述，基于DEM提取地性線是項(xiàng)復(fù)雜的工作，目前提取算法相對成熟，但不同算法各有利弊，針對實(shí)際工作的不同地域地形和收集到的DEM數(shù)據(jù)質(zhì)量，應(yīng)選擇最佳的算法進(jìn)行提取。