王 璐，艾廷華

武漢大學(xué)資源與環(huán)境科學(xué)學(xué)院，湖北 武漢 430079

針對該問題，本文探討了六邊形DEM構(gòu)建以及谷地線自動提取在六邊形DEM格網(wǎng)結(jié)構(gòu)上的實現(xiàn)，并以等高線數(shù)據(jù)分別建立四邊形DEM和六邊形DEM，采用多分辨率分析方法，以谷地線數(shù)量、長度和形狀3個特征指標(biāo)對2種DEM所提取的結(jié)果進(jìn)行討論，以此評價DEM網(wǎng)格幾何形態(tài)對于谷地線提取的影響。

1 六邊形格網(wǎng)DEM建立

1.1 六邊形格網(wǎng)結(jié)構(gòu)

格網(wǎng)結(jié)構(gòu)是DEM表達(dá)真實地理空間的數(shù)字格網(wǎng)空間基礎(chǔ)，為了滿足DEM應(yīng)用分析，六邊形格網(wǎng)結(jié)構(gòu)的建立應(yīng)包括：①建立適合計算機存儲與處理，易于與矢量坐標(biāo)實現(xiàn)對接的格網(wǎng)坐標(biāo)系；②建立格網(wǎng)之間的鄰域關(guān)系。

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

圖1 六邊形格網(wǎng)結(jié)構(gòu)Fig.1 Hexagon gridded structure

(1)

式中，X(i,j)和Y(i,j)為格網(wǎng)坐標(biāo)為(i,j)網(wǎng)格中心點的地理坐標(biāo)；x0為格網(wǎng)坐標(biāo)原點所對應(yīng)的地理橫坐標(biāo)；y0為格網(wǎng)坐標(biāo)原點所對應(yīng)的地理縱坐標(biāo)。

6鄰域關(guān)系是六邊形格網(wǎng)結(jié)構(gòu)應(yīng)用于DEM分析中鄰域處理的基本關(guān)系。然而，行之間的錯位使奇數(shù)行格網(wǎng)和偶數(shù)行格網(wǎng)的6鄰域格網(wǎng)坐標(biāo)編碼規(guī)律存在差異，如圖1(b)標(biāo)紅鄰域，偶數(shù)行格網(wǎng)中該鄰域的列號均比奇數(shù)行格網(wǎng)大1，以此編碼規(guī)律便可構(gòu)建格網(wǎng)之間的鄰域關(guān)系。

1.2 規(guī)則格網(wǎng)DEM構(gòu)建

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

(2) 根據(jù)格網(wǎng)坐標(biāo)與地理坐標(biāo)之間的對應(yīng)關(guān)系，獲取格網(wǎng)結(jié)構(gòu)中每個格網(wǎng)單元中心點對應(yīng)的地理坐標(biāo)，以便后續(xù)格網(wǎng)高程值的內(nèi)插計算。

(3) 利用等高線構(gòu)建TIN數(shù)據(jù)結(jié)構(gòu)，建立六邊形格網(wǎng)中心點與TIN三角形的對應(yīng)關(guān)系(圖2)，對于每個格網(wǎng)中心點，獲取其所在TIN結(jié)構(gòu)中的三角形，基于三角形的線性內(nèi)插法計算格網(wǎng)中心點的高程值。

圖2 六邊形格網(wǎng)中心點與TIN三角形對應(yīng)Fig.2 Correspondence between central points of hexagons and triangles of TIN

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

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

2.1 填洼處理

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

該算法應(yīng)用到六邊形DEM格網(wǎng)結(jié)構(gòu)中原理相同，而最大的不同之處在于由出水口向內(nèi)部搜索時對于格網(wǎng)單元由8方向到6方向的鄰域處理變化，具體實現(xiàn)步驟為：

(1) 將六邊形DEM數(shù)據(jù)的邊緣格網(wǎng)高程值設(shè)置為溢出高程值，并壓入最小優(yōu)先隊列，將相應(yīng)格網(wǎng)標(biāo)記為“1”。

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

Hn=max{Ec,En}

(2)

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

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

(4) 重復(fù)步驟(2)、(3)，直至最小優(yōu)先隊列為空，填平結(jié)束。

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

2.2 D6算法

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

(1) 格網(wǎng)水流方向編碼。根據(jù)六邊形格網(wǎng)結(jié)構(gòu)中鄰域關(guān)系，規(guī)定水流從中心格網(wǎng)流向正東方向格網(wǎng)為初始方向，逆時針對該6個鄰域單元流向編碼為1，2，…，6，該編碼也為鄰域格網(wǎng)順序編號。

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

Gn=zc-zn

(3)

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

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

圖3 六邊形DEM格網(wǎng)結(jié)構(gòu)下W&L填洼算法流程Fig.3 W&L depression filling algorithm based on hexagon gridded structure

2.3 平地區(qū)域水流方向計算

經(jīng)過填洼處理，D6算法流向計算結(jié)果分為3種情況：①僅有一個具有最大坡降值的鄰域格網(wǎng)，如圖4(a)，該情況下水流方向值唯一；②多個具有最大坡降值大于0的鄰域格網(wǎng)，如圖4(b)，該情況下可按編碼順序?qū)⑺鞣较蛑赶虻谝粋€鄰域格網(wǎng)；③多個具有最大坡降值等于0的鄰域格網(wǎng)，如圖4(c)，該情況下為平地格網(wǎng)，若按照情況②進(jìn)行處理，易出現(xiàn)“流向互指”現(xiàn)象，影響后續(xù)計算，需對平地格網(wǎng)進(jìn)行特殊計算。

圖4 D6算法Fig.4 D6 algorithm

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

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

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

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

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

步驟3，重復(fù)步驟2，直至待處理數(shù)組為空，平地區(qū)域水流方向計算結(jié)束。

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

圖5 平地格網(wǎng)水流計算Fig.5 Flow direction calculation for flat area grids

基于六邊形格網(wǎng)的鄰域處理，結(jié)合W&L填洼算法、D6算法以及平地區(qū)域流向的計算為基于六邊形DEM谷地線的提取提供明確的水流方向，如圖6，后續(xù)水流累計量計算依賴于流向結(jié)果完成，提取累積量大于2的格網(wǎng)單元，并按照對應(yīng)格網(wǎng)的流向線作拓?fù)溥B通組織，完成谷地線的追蹤提取。

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

3 試驗與分析

3.1 規(guī)則格網(wǎng)DEM建立結(jié)果

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

表1 兩種DEM高程中誤差計算結(jié)果

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

圖7 四邊形DEM構(gòu)建結(jié)果Fig.7 The results of square gridded DEM generation

圖8 六邊形DEM構(gòu)建結(jié)果Fig.8 The results of hexagon gridded DEM generation

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

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

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

3.2 谷地線提取結(jié)果對比分析

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

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

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

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

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

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

(5)

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

圖12 誤差帶Fig.12 Error band

4 結(jié) 論

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

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