正六邊形規則格網表達的DEM谷地線提取

王 璐，艾廷華

武漢大學資源與環境科學學院，湖北 武漢 430079

針對該問題，本文探討了六邊形DEM構建以及谷地線自動提取在六邊形DEM格網結構上的實現，并以等高線數據分別建立四邊形DEM和六邊形DEM，采用多分辨率分析方法，以谷地線數量、長度和形狀3個特征指標對2種DEM所提取的結果進行討論，以此評價DEM網格幾何形態對于谷地線提取的影響。

1 六邊形格網DEM建立

1.1 六邊形格網結構

格網結構是DEM表達真實地理空間的數字格網空間基礎，為了滿足DEM應用分析，六邊形格網結構的建立應包括：①建立適合計算機存儲與處理，易于與矢量坐標實現對接的格網坐標系；②建立格網之間的鄰域關系。

正六邊形格網坐標系包括正交行列坐標系、極坐標系、三斜軸坐標系以及Gosper曲線坐標系[15-16]。其中正交行列坐標系是標記六邊形格網單元最簡單的方式，其思想與四邊形格網結構類似，即用行列號直接對格網單元進行編碼，該坐標系與笛卡兒坐標系兼容，便于計算機存儲與組織，易于實現從格網坐標到地理坐標的轉換。本文基于“奇數行左偏移”的正交行列坐標系構建六邊形格網結構，如圖1所示，以(i,j)為對格網坐標進行編碼，其中i和j分別為格網單元所在的列、行數，(0,0)格網中心點為坐標系原點，該結構中行與行之間呈“之”字形錯位。

圖1 六邊形格網結構Fig.1 Hexagon gridded structure

(1)

式中，X(i,j)和Y(i,j)為格網坐標為(i,j)網格中心點的地理坐標；x0為格網坐標原點所對應的地理橫坐標；y0為格網坐標原點所對應的地理縱坐標。

6鄰域關系是六邊形格網結構應用于DEM分析中鄰域處理的基本關系。然而，行之間的錯位使奇數行格網和偶數行格網的6鄰域格網坐標編碼規律存在差異，如圖1(b)標紅鄰域，偶數行格網中該鄰域的列號均比奇數行格網大1，以此編碼規律便可構建格網之間的鄰域關系。

1.2 規則格網DEM構建

為保證六邊形DEM和四邊形DEM在建立過程中誤差來源的一致性，本文利用等高線構建規則格網DEM[20]。

(2) 根據格網坐標與地理坐標之間的對應關系，獲取格網結構中每個格網單元中心點對應的地理坐標，以便后續格網高程值的內插計算。

(3) 利用等高線構建TIN數據結構，建立六邊形格網中心點與TIN三角形的對應關系(圖2)，對于每個格網中心點，獲取其所在TIN結構中的三角形，基于三角形的線性內插法計算格網中心點的高程值。

圖2 六邊形格網中心點與TIN三角形對應Fig.2 Correspondence between central points of hexagons and triangles of TIN

2 基于六邊形DEM的谷地線提取

谷地線提取的原理[7]源于地表徑流的自然特性，即水流沿斜坡最陡方向，通過比較鄰域格網的高程值計算各格網點水流方向，基于水流方向計算匯合到每一個格網的地表徑流量，通過設定流量閾值，提取大于該閾值的格網即可獲得相應的谷地線起始位置以及谷地線網絡。簡而言之，谷地線自動提取過程包括：填洼處理、水流方向計算、水流累積量計算及流量閾值確定，最后根據流向追蹤提取谷地線。其中填洼處理、水流方向計算及平地區域水流方向確定決定了后續步驟的計算，為谷地線提取的關鍵問題，為此，國內外提出了眾多算法[21-23]，然而算法的實現均基于四邊形DEM格網結構以及鄰域關系，未對六邊形DEM格網結構下算法實現的可行性進行討論。因此，本文基于W&L算法、單流向算法分別討論在六邊形格網DEM結構中填洼處理、流向計算的實現，并為平地區域設計流向計算方法。

2.1 填洼處理

W&L算法是高效率填洼算法之一[24]。通過模擬洪水淹沒自然地形的過程，即水平面從最低點逐漸抬升至淹沒整個地表，柵格單元從外圍向內部逐個被淹沒，該算法以DEM的邊緣最低點作為出水口點，通過與鄰域格網點高程值進行比較，以確定每個格網單元流向出水口的最小高程值為溢水高程值，并將格網單元高程值提升至溢水高程值，其實質就是從潛在出口向內部進行逆向鄰域搜索完成整個填洼過程。其實現過程是在溢水高程值的基礎上，利用最小代價搜索，以優先隊列數據結構完成填洼。

該算法應用到六邊形DEM格網結構中原理相同，而最大的不同之處在于由出水口向內部搜索時對于格網單元由8方向到6方向的鄰域處理變化，具體實現步驟為：

(1) 將六邊形DEM數據的邊緣格網高程值設置為溢出高程值，并壓入最小優先隊列，將相應格網標記為“1”。

(2) 獲取最小優先隊列中隊首元素以及對應的格網單元c，根據最小鄰域單元查找6鄰域格網中所有未被標記的格網，并利用下式計算鄰域格網的溢出高程值

Hn=max{Ec,En}

(2)

式中，Hn為第n個鄰域格網的溢出高程值；Ec和En分別為格網單元c以及第n個鄰域格網的原始高程值。

(3) 從最小優先隊列中刪除格網單元c對應的溢出高程值，將步驟(2)中未被標記的鄰域格網溢出高程值壓入最小優先隊列，并將其標注為已處理。

(4) 重復步驟(2)、(3)，直至最小優先隊列為空，填平結束。

圖3為算法示意圖，第2次迭代中，C為隊首元素所對應的格網單元，通過查找可知其鄰域單元中F和G未被標記，根據計算公式可得F和G的溢出高程均為40，將F和G壓入最小優先隊列，刪除C，即完成了F和G格網的填洼。同理，格網K和J分別在第3、4次迭代中被填平，直至第16次迭代，整個區域完成填洼過程。根據上述算法實現步驟，W&L算法原理應用到六邊形格網結構中的鄰域關系，也可為六邊形DEM實現填洼處理。

2.2 D6算法

類比于四邊形DEM中應用最為廣泛的單流向算法——D8算法[7]，本文提出D6算法以計算六邊形DEM中格網的流向。該算法基于最大坡降法，可定義為：對于任何一個六邊形格網c，在其相鄰六個格網點中選擇滿足高程落差值最大的鄰域格網點作為其流向(即水流的流出方向)。其具體計算步驟包括：

(1) 格網水流方向編碼。根據六邊形格網結構中鄰域關系，規定水流從中心格網流向正東方向格網為初始方向，逆時針對該6個鄰域單元流向編碼為1，2，…，6，該編碼也為鄰域格網順序編號。

(2) 鄰域格網坡降計算。根據式(3)，分別計算中心格網高程值與其6個鄰域格網的高程值差值分別作為其鄰域格網坡降值

Gn=zc-zn

(3)

式中，Gn為中心格網在編號為n的鄰域格網方向上的坡降值；zc和zn分別為中心格網和編號為n的鄰域格網的高程值。

(3) 水流方向確定。根據(2)中的計算結果，令中心格網水流流向具有最大坡降值的鄰域格網。

圖3 六邊形DEM格網結構下W&L填洼算法流程Fig.3 W&L depression filling algorithm based on hexagon gridded structure

2.3 平地區域水流方向計算

經過填洼處理，D6算法流向計算結果分為3種情況：①僅有一個具有最大坡降值的鄰域格網，如圖4(a)，該情況下水流方向值唯一；②多個具有最大坡降值大于0的鄰域格網，如圖4(b)，該情況下可按編碼順序將水流方向指向第一個鄰域格網；③多個具有最大坡降值等于0的鄰域格網，如圖4(c)，該情況下為平地格網，若按照情況②進行處理，易出現“流向互指”現象，影響后續計算，需對平地格網進行特殊計算。

圖4 D6算法Fig.4 D6 algorithm

對于無流向值平地格網，通過六鄰域搜索尋找流向出水口的路徑，并由該路徑逆序計算平地格網流向值，具體方法為：

(1) 基于D6算法確定非平地格網單元流向值，將平地格網單元存入待處理數組中。

(2) 對待處理數組中每個平地格網作如下處理：

步驟1，按順序查找6鄰域中高程值與中心格網相等的鄰域格網，若存在已有流向值且流向不指向中心格網的鄰域格網，則將中心格網的流向指向該鄰域，并刪除該格網。

步驟2，若不存在步驟1中的鄰域格網，則依次以無流向值的鄰域格網為中心格網，執行步驟1。

步驟3，重復步驟2，直至待處理數組為空，平地區域水流方向計算結束。

圖5中的平地區域是由2.1節中W&L算法填洼處理后所形成，其中格網C、F、G、J和K為平地格網。以格網C流向計算為例，其鄰域格網中不存在已有流向值且流向不指向中心格網的鄰域格網，因此按順序對其無流向值的鄰域格網F和G遍歷，第1次F鄰域搜索中未找到流向出水口的路徑，第2次G鄰域搜索中查找到格網L為出水口P的相鄰格網，則路徑為C→G→L→P，逆序對G和C流向進行計算，結果均為6。同理，計算格網F、J和K流向分別為1、2和1。由計算結果可知，該方法可有效避免平地區域流向結果“互指”現象。

圖5 平地格網水流計算Fig.5 Flow direction calculation for flat area grids

基于六邊形格網的鄰域處理，結合W&L填洼算法、D6算法以及平地區域流向的計算為基于六邊形DEM谷地線的提取提供明確的水流方向，如圖6，后續水流累計量計算依賴于流向結果完成，提取累積量大于2的格網單元，并按照對應格網的流向線作拓撲連通組織，完成谷地線的追蹤提取。

圖6 谷地線提取Fig.6 Extraction of valley lines

3 試驗與分析

3.1 規則格網DEM建立結果

本文試驗區域為四川省德陽市西北部龍門山脈中段，該區域地勢西北高東南低，谷地發育明顯。試驗數據為國家測繪部門標準化生產的1∶5萬DLG等高線數據，其等高距為10 m，最高高程為4400 m，最低高程為1000 m，以此數據分別構建九種分辨率的四邊形DEM和六邊形DEM數據，相同分辨率下四邊形DEM和六邊形DEM水平方向格網寬度比值約為0.93，見表1。

表1 兩種DEM高程中誤差計算結果

如圖7和圖8，高分辨率下，基于六邊形格網結構和四邊形格網結構所構建的DEM對地形可視化效果差異不大，而在低分辨率下，不同高程值格網之間的界線以及區域邊界在六邊形DEM中表現更為明顯、平滑，地形變化更直觀。

圖7 四邊形DEM構建結果Fig.7 The results of square gridded DEM generation

圖8 六邊形DEM構建結果Fig.8 The results of hexagon gridded DEM generation

為保證后續谷地線提取結果的可信性，本文采用檢查點法[25]，隨機選取研究區域內和邊緣的42個地貌特征點作為檢測點，以格網點高程中誤差對所建立的六邊形DEM和四邊形DEM數據精度分別進行評定，其公式為

式中，RMSE為DEM高程中誤差；Zk(k=1,2,…,n)為檢測點高程值；Rk為DEM內插所得對應檢測點高程值。

由表1結果可知，本文基于等高線所建立的2種DEM在分辨率1的高程中誤差分別為9.33 m和9.37 m，分辨率2的高程中誤差為11.29 m和10.55 m。根據我國1∶5萬DEM精度標準中對于高山地區域25 m格網間隔DEM的中誤差限差為19 m的要求[25]，以本文方法所建立的四邊形格網DEM在格網間隔為25 m時，精度在9.33 m至11.29 m之間，與該分辨率對應的六邊形DEM的格網寬度約為27 m，精度在9.37 m至10.55 m之間，均符合高程中誤差限差19 m的要求。此外，隨著分辨率降低，所構建的六邊形DEM精度損失程度整體上小于四邊形DEM，尤其在分辨率8和9時，其精度損失程度遠低于四邊形DEM。

3.2 谷地線提取結果對比分析

根據3.1節中構建的四邊形DEM和六邊形DEM，按照谷地線提取步驟，基于D8算法和D6算法，以25為累積流量閾值分別提取不同分辨率下的谷地線網絡。如表2所示，相同分辨率下，六邊形DEM所提取的谷地線長度均大于四邊形DEM提取結果，但谷地線數量差異并不大。以1∶5萬DLG水系數據作為參考水系，將提取結果分別與參考水系進行對比可知，分辨率1時，兩種DEM所提取的谷地線主干大致相同，但在地形起伏較小區域均易出現“平行谷地線”，如圖9中綠框區域，但兩者鄰域之間的角度不同使該區域谷地線的流向存在差異，六邊形DEM谷地線流向與正東方向所成夾角多呈30°，而四邊形DEM則呈45°或90°。

表2 兩種DEM所提取的谷地線特征比較

隨著分辨率減小，兩種DEM下所提取的谷地線長度和數量均隨之減少，分支丟失，形狀逐漸被簡化，但兩者提取結果差異變大。低分辨率，六邊形DEM在保持谷地線形狀彎曲特征上效果更好，四邊形DEM下提取的谷地線形狀簡化程度明顯，尤其是對于縱向支流的提取，如圖9中藍框區域。為了更直觀的對比該區域中不同分辨率下兩者提取結果的簡化程度，將提取結果進行放大顯示，由圖10和圖11可知，分辨率2，兩者結果形狀簡化均不明顯，數量和長度簡化較為明顯，四邊形DEM所提取的谷地線在分辨率5時趨于平直，分辨率9，被簡化為直線，反觀六邊形DEM提取結果仍保持一定的彎曲度。此外，分辨率9時，基于四邊形DEM的谷地線提取結果在圖9(d)中紅色區域中出現斷裂的情況，而六邊形DEM提取結果仍保持水流線之間相互連接關系。

圖10 D8算法谷地線提取結果局部放大圖Fig.10 Enlarged results of extracted valley lines with D8 algorithm

圖11 D6算法谷地線提取結果局部放大圖Fig.11 Enlarged results of extracted valley lines with D6 algorithm

為了定量化對比2種DEM所提取谷地線形狀特征，本文以誤差帶寬度表示各分辨率下谷地線提取結果偏離參考水系的程度

(5)

式中，Di為分辨率i下誤差帶寬度；Δs為套合面面積；Li為分辨率i下參考數據套合線長度。Di值越小， 所提取谷地線形狀特征與參考水系吻合程度越高。由表2誤差帶寬度結果可知，隨著分辨率減小，四邊形DEM和六邊形DEM所提取的谷地線誤差帶寬度呈現出先減小后增大的趨勢，均在分辨率3時出現最小值，分別為29.88 m和29.99 m，即該分辨率下兩者提取結果與參考水系最為吻合。分辨率4—9，兩種DEM所提取的谷地線與參考水系形成的誤差帶寬度隨分辨率減小而增大，吻合程度降低，但誤差帶寬度在六邊形DEM中均低于四邊形DEM，分別低了6.75%、2.76%、5.64%、6.33%、6.13%及2.94%，這表明，中、低分辨率下六邊形DEM提取結果與參考水系吻合程度較高，在保持形狀特征方面優于四邊形DEM，較大程度的保持了彎曲特征。

圖12 誤差帶Fig.12 Error band

4 結 論

本文從DEM結構角度，采用多分辨率分析方法，探討格網幾何形態對于谷地線提取的影響。鑒于六邊形具有鄰域一致性、各向同性、緊湊性、采樣率高等優點，本文以六邊形為基本格網單元構建DEM，并基于六鄰域處理單元，對六邊形DEM谷地線提取過程中填洼、流向計算以及平地區域流向確定三個關鍵步驟進行討論實現，從而實現谷地線在六邊形DEM中的提取。試驗結果表明：在DEM構建方面，六邊形格網結構可保持更高的數據精度；在谷地線提取方面，六邊形格網結構的優勢表現為對于谷地線形狀的保持，隨著分辨率減小，六邊形DEM所提取的谷地線與實測數據吻合程度高，彎曲特征繼承性強。因此，在相同數據存儲量條件下，六邊形DEM數據精度更高，且其提取的谷地線網絡的形狀特征更為精細。

然而，本文集中于討論四邊形與六邊形在谷地線提取的差異，未考慮到提取算法的優化，仍采用單流向D6算法，其提取結果仍無法避免“平行谷地線”，使其與真實谷地線不符，對此，后續研究還需結合多流向算法或利用真實谷地線網絡對DEM進行修正從而改善“平行谷地線”現象。此外，從應用角度，進一步的研究包括將六邊形DEM谷地線應用在后續流域劃分、匯流面積計算等水文參數的提取中。