先交代一下背景這套流程是我在做一個流域地形分析項(xiàng)目時反復(fù)打磨出來的。ArcGIS空間分析模塊里的水文分析工具絕大多數(shù)人只用它提取河網(wǎng)、劃分子流域很少有人會想到把DEM反過來用一遍直接提取山谷線和山脊線。這篇文章就把整套方法和盤托出從核心思路到每一步操作細(xì)節(jié)再到我踩過的一些坑一次說清楚。適合正在做地形分析、從事GIS數(shù)據(jù)生產(chǎn)、或者剛接觸空間分析但對算法原理還不太明白的朋友閱讀按文章順序操作一遍你也能在ArcGIS里穩(wěn)定提取出山谷線和山脊線。1. 先把核心思路理清楚山谷線是“流水路”山脊線是“分水墻”1.1 山谷線和山脊線在地貌上到底“長什么樣”在動手之前先建立直觀認(rèn)識。山谷線也叫溝谷線是降雨后地表水最容易匯集的連線。你在地形圖上看等高線等高線向高處凸出的那串鞍部連線基本就是山谷線放到真實(shí)地貌里就是山溝、溪流所在的位置。山脊線則剛好相反它是分水嶺雨水落在山脊線上會往兩側(cè)分開流地形上表現(xiàn)為等高線向低處凸出的連線也就是山梁的頂部走向。這兩類線是地形分析里最基本的地形骨架線。做水文分析要劃分子流域分水嶺山脊線是天然邊界做道路選線要盡量沿山脊或緩坡布線需要知道山脊線的精確位置做植被分布調(diào)查、地質(zhì)災(zāi)害評估也經(jīng)常要參考溝谷和山梁的走向?？蓡栴}在于人工勾繪山脊山谷線非常費(fèi)勁尤其研究區(qū)一大等高線密密麻麻眼睛看花不說不同人勾出來的結(jié)果還不一樣標(biāo)準(zhǔn)很難統(tǒng)一。如果只看一句水文上的定義那就是山谷線是匯水線山脊線是分水線。這句話是整個方法的底層依據(jù)。1.2 地形翻轉(zhuǎn)用一次反向DEM把山脊變成“假河谷”ArcGIS的水文分析工具本身是為提取河網(wǎng)設(shè)計的。水流在DEM上從高往低走山谷線天然滿足“水流匯聚”這個條件所以直接用Flow Direction計算流向、再用Flow Accumulation計算匯流累積量就能把山谷線找出來。問題出在山脊線上。山脊線是水流的源頭、是分水的邊界水流不會沿著山脊走反而會從山脊向兩側(cè)散開。你用同樣的工具去分析會發(fā)現(xiàn)山脊線上的匯流累積值非常小根本不可能作為“河網(wǎng)”被提取出來。解決辦法其實(shí)是一招思維反轉(zhuǎn)把DEM的每個像元乘以-1山變谷、谷變山。原本的山脊線在反向DEM上變成了“山谷”水流會在這里匯聚于是一套完全相同的流向、流量分析流程跑下來提取出的“反向河網(wǎng)”就是山脊線。這個思路并不復(fù)雜但理解以后你會發(fā)現(xiàn)它的可擴(kuò)展性很強(qiáng)。提取溝沿線、提取鞍部點(diǎn)、分析地形起伏特征本質(zhì)上都是在跟“地形翻轉(zhuǎn)”做文章。把底層邏輯想明白后面調(diào)參數(shù)時就不會慌。2. 開始前的準(zhǔn)備工作不是隨便拿個DEM就能跑2.1 DEM數(shù)據(jù)的要求與常見坑理論上來講任何DEM都能跑這套流程但結(jié)果的可靠程度取決于DEM質(zhì)量和分辨率。分辨率一般用10米到30米的DEM比較合適。分辨率太高比如1米或2米地形細(xì)節(jié)太碎提取出的線會非常細(xì)碎后期要花大量時間修整分辨率太低比如90米或250米山谷和山脊的輪廓被壓縮很多次級溝谷根本顯示不出來。數(shù)據(jù)源國內(nèi)做研究常用地理空間數(shù)據(jù)云下載的ASTER GDEM、GDEMV3或者ALOS 12.5米數(shù)據(jù)做工程應(yīng)用時我一般用項(xiàng)目自己生成的5米或2米DEM。投影坐標(biāo)系這是很多人忽略的坑。水文分析里的流向計算基于平面距離如果數(shù)據(jù)是經(jīng)緯度坐標(biāo)系每個像元在不同緯度上的實(shí)際距離不一樣流向判斷會出問題。實(shí)際操作時先把DEM投影到合適的坐標(biāo)系比如UTM或者CGCS2000高斯投影再做水文分析結(jié)果才可靠。NoData區(qū)域DEM邊緣或者水域位置如果存在NoData需要提前處理。流向計算遇到NoData會把像元當(dāng)成阻礙導(dǎo)致河網(wǎng)在邊界處亂跑甚至直接斷裂。另外提醒一句拿到DEM先檢查屬性表里的像元大小和范圍。ArcGIS的柵格分析對一致性非常敏感后面所有工具是否能順利跑通很大程度上取決于原始數(shù)據(jù)干不干凈。2.2 Spatial Analyst擴(kuò)展模塊與工具集位置水文分析這組工具位于ArcToolbox里的Spatial Analyst Tools下的Hydrology工具集。使用前要確保ArcMap或ArcGIS Pro里已經(jīng)啟用Spatial Analyst擴(kuò)展模塊。在ArcMap中路徑是“自定義”菜單里的“擴(kuò)展模塊”勾選Spatial Analyst。在ArcGIS Pro中從“設(shè)置”進(jìn)入“許可”頁面啟用Spatial Analyst。沒有這個許可很多柵格分析工具會直接灰掉或者報錯。企業(yè)版許可一般自帶但個人使用時要特別注意別等到跑到一半才發(fā)現(xiàn)工具不可用。順便說明一下ArcGIS 10.x中這些工具的名稱和ArcGIS Pro 3.x基本一致區(qū)別主要在界面位置。ArcMap里在ArcToolbox中操作Pro里在“分析”選項(xiàng)卡下的“工具”面板中查找。下面講的流程在兩個平臺通用差別只是面板叫法不同。3. 正向水文分析提取山谷線從填洼到矢量化的完整流程3.1 填洼先把DEM里的“坑”處理明白水文分析的第一步是填洼。DEM里因?yàn)閿?shù)據(jù)噪聲或真實(shí)地形原因會存在一些低洼區(qū)域。如果不填水流流到洼地就停住了后續(xù)流向計算會斷掉河網(wǎng)自然提不出來。填洼工具在Hydrology工具集下名稱就叫Fill。ArcMap里直接輸入DEM就能跑默認(rèn)參數(shù)會把所有洼地都填平。這里有一個容易被忽略的Z limit參數(shù)如果留空表示所有洼地都會被填掉。但實(shí)際地形中確實(shí)存在真實(shí)洼地比如喀斯特地區(qū)、采石坑、水庫這些地方一味填平會把真實(shí)地形破壞掉導(dǎo)致提取的河網(wǎng)出現(xiàn)“穿過山體”的假象。我的習(xí)慣做法是先用Fill做一次全填然后用柵格計算器求原始DEM和填洼后DEM的差值差值大于一定深度的區(qū)域就是重點(diǎn)關(guān)注對象。如果研究區(qū)不大可以先把這些差值區(qū)域目視檢查一遍再決定要不要限制Z limit。比如某個區(qū)域差值達(dá)到50米以上那很可能是真實(shí)洼地或人工地物填掉它會影響后續(xù)分析。這一步還有個小細(xì)節(jié)如果DEM邊緣有NoData建議先用“按掩膜提取”把研究區(qū)邊界裁好或者使用“消除NoData”處理掉邊緣無效值再進(jìn)入填洼流程。否則填洼工具在NoData邊緣可能會產(chǎn)生不可預(yù)期的結(jié)果。3.2 流向計算D8算法給每個像元指路填洼完成后用Flow Direction工具計算流向。ArcGIS默認(rèn)使用D8單流向算法它用一個3x3的窗口計算中心像元到周圍8個鄰域像元的坡度選坡度最大的方向作為該像元的水流去向。參數(shù)上重點(diǎn)說一下“流向類型”默認(rèn)是D8適合河網(wǎng)提取。至于MF-D8多流向方向在平原區(qū)效果更好但后面做流量累積時算法需要配套調(diào)整不建議在基礎(chǔ)流程里使用。D-Infinity連續(xù)流向多用于徑流擴(kuò)散研究也不常用。另一個選項(xiàng)Force all edge cells to flow outward意思是強(qiáng)制邊緣像元向外流。默認(rèn)不勾選因?yàn)槿绻鸇EM邊緣就是研究區(qū)邊界強(qiáng)行讓邊緣水流向外會生成排水方向異常的柵格。但如果你處理的是分塊數(shù)據(jù)邊緣確實(shí)還有地形延續(xù)可以考慮勾選。流向柵格生成后我習(xí)慣順手檢查一下它的屬性表。注意流向柵格不是簡單記錄方向數(shù)值而是用2的冪次方編碼表示8個方向分別是1、2、4、8、16、32、64、128。在可視化時你可以用唯一值渲染來看效果如果出現(xiàn)大片方向值異常大概率是前面填洼或DEM處理環(huán)節(jié)出了問題。3.3 流量累積與河網(wǎng)閾值數(shù)字越大線越少Flow Accumulation工具的輸入是流向柵格輸出結(jié)果表示每個像元上游累積了多少個像元的匯水量。比如某個像元的結(jié)果是1000意味著它上游有1000個像元的水最終匯到這個位置。在沒有雨量權(quán)重的情況下可以把這個值理解為“上游匯水面積”。提取河網(wǎng)時關(guān)鍵環(huán)節(jié)是設(shè)置一個閾值把流量累積超過閾值的像元定義為中心線。那么閾值怎么選經(jīng)驗(yàn)上先用流量累積最大值的1%作為起點(diǎn)。例如流量累積柵格最大值是50000就先試500。然后疊加到影像或等高線上看效果。閾值設(shè)得越低提取的溝谷越密集閾值越高只保留主溝谷。沒有一個絕對正確的值取決于你要表達(dá)的是“完整山溝系統(tǒng)”還是“主要河谷”。實(shí)際工作中我通常先用一個比較低的閾值跑一版疊加到底圖上看。如果發(fā)現(xiàn)很多“溝”其實(shí)只是斜坡上的細(xì)碎沖溝就調(diào)高閾值再跑。流量累積結(jié)果往往呈長尾分布主河谷區(qū)域值很大支溝的值相對小。如果你發(fā)現(xiàn)不管怎么調(diào)閾值主線附近的支溝總是亂糟糟的那可能是原始DEM里含有大量噪聲比如林區(qū)雷達(dá)數(shù)據(jù)生成的DEM有樹冠噪聲這種情況下需要先對DEM做平滑濾波處理。3.4 柵格河網(wǎng)轉(zhuǎn)矢量線Stream Order與Stream to Feature得到閾值后的河網(wǎng)柵格后我用柵格計算器寫Con(flowacc 閾值, 1)把河網(wǎng)像元賦為1背景設(shè)為NoData。這一步也可以直接使用SetNullflowacc 閾值, flowacc實(shí)現(xiàn)。然后有兩個選擇。一種是把河網(wǎng)柵格做Stream Order分級再進(jìn)行矢量化。Stream Order有兩種常用分級方法Strahler和Shreve。Strahler分級是主流它把最上游支流定為1級兩條同級支流交匯后級別加1。這樣你可以把不同級別的河谷線分別存儲、分別顯示線狀符號按級別區(qū)分粗細(xì)圖面效果非常專業(yè)。另一種是直接用Stream to Feature工具。這個工具在實(shí)際調(diào)用時會自動執(zhí)行Stream Order但前提是你的河網(wǎng)柵格像元值干凈。我的建議是先把河網(wǎng)柵格跑一遍Stream Order再用Stream to Feature輸入分級結(jié)果這樣輸出線的屬性表里會帶有Strahler等級字段gridcode方便后續(xù)按等級篩選和符號化。勾選Simplify選項(xiàng)可以去掉柵格轉(zhuǎn)矢量時產(chǎn)生的鋸齒但簡化過度會犧牲真實(shí)位置精度一般保留默認(rèn)設(shè)置后面再用Smooth Line做一次平滑。矢量化出來的線就是山谷線。4. 反向地形提取山脊線關(guān)鍵就在“DEM乘以-1”4.1 反地形DEM生成注意數(shù)據(jù)類型別被截斷反向DEM其實(shí)就是地形翻轉(zhuǎn)。在ArcGIS的柵格計算器里直接寫DEM * -1或者使用“數(shù)學(xué)分析”工具里的“減”0 - DEM。這里有一個容易被忽視的問題反向后DEM的取值范圍從原來的正值變成負(fù)值這沒關(guān)系流向計算只關(guān)心相對高差。但如果你后續(xù)還有其他分析要用到高程排序負(fù)值會影響邏輯判斷。還有一個數(shù)據(jù)格式陷阱。如果DEM是整型柵格直接乘以-1的結(jié)果可能仍然是整型像元值被截斷精度會損失。所以我在柵格計算器里通常先加一步轉(zhuǎn)浮點(diǎn)用Float函數(shù)處理一下確保精度不丟。這個細(xì)節(jié)在DEM分辨率較高時影響不大但在地勢平緩區(qū)域微小的高程差異可能決定一條山脊線的走向。做到這一步原本高處的山脊變成了洼地原本深谷變成了高地整個地形被翻轉(zhuǎn)過來。接下來對反向DEM再做一次完整的填洼、流向、流量累積流程。4.2 反向水文分析填洼、流向、累積全流程細(xì)節(jié)操作順序和正向分析完全一樣Fill反向DEM→ Flow Direction → Flow Accumulation → 閾值 → 柵格轉(zhuǎn)矢量。但有幾個細(xì)節(jié)值得單獨(dú)拆開說。首先反向填洼非常關(guān)鍵。這個流程的本質(zhì)是找反向DEM中的“匯水線”如果填洼過度很多真正反映山脊形態(tài)的細(xì)節(jié)洼地被填掉提取出的山脊線會明顯偏離真實(shí)山脊往往偏向山脊的任一側(cè)。解決方法是填洼時限制Z limit或者比較填洼前后的差異選擇性地保留部分負(fù)地形。其次反向河網(wǎng)的閾值通常比正向河網(wǎng)難調(diào)。因?yàn)樯郊咕€地形比山谷線地形更破碎同樣的閾值下支脊數(shù)量會多很多。實(shí)際工作中我傾向于“主脊優(yōu)先”策略先把閾值調(diào)高提取主要山脊線再逐步降低閾值看次級分支而不是一上來就追求把所有細(xì)碎山脊一次提取完。流向計算完成后我建議把反向流量累積圖疊加到山體陰影或者等高線上目視驗(yàn)證一次。真正的山脊線在反向流量累積圖上應(yīng)該是連續(xù)的亮色線位置基本與等高線轉(zhuǎn)折、山體陰影的分界線吻合。如果發(fā)現(xiàn)提取的線走在山坡中間多半是填洼深度過大或者原始DEM存在偽洼地。4.3 山脊線的后處理平滑、刪短小偽線反向河網(wǎng)柵格轉(zhuǎn)換矢量的操作跟正向一樣用Stream Order Stream to Feature處理得到的線就是山脊線。由于反向DEM在地形翻轉(zhuǎn)后相鄰像元間的微小起伏會被放大提取的矢量線通常比山谷線更曲折、毛刺更多。所以矢量化之后我一般會再執(zhí)行一次Smooth Line平滑工具。平滑使用PAEK算法容差可以設(shè)為DEM像元尺寸的2到5倍。比如30米DEM平滑容差設(shè)為60到150米。容差太大會把山脊線的關(guān)鍵拐點(diǎn)削平容差太小又沒有實(shí)際效果。另一個優(yōu)化是在柵格階段做。Flow Accumulation工具支持輸入權(quán)重柵格如果你手里有坡度數(shù)據(jù)可以把坡度作為權(quán)重讓陡坡區(qū)域更容易累積出高流量從而突出陡峻山脊。不過這個方法不一定適合所有地形如果只是想快速出圖可以先不加權(quán)重。在矢量階段我還會根據(jù)線長篩選。提取出的山脊線中有些很短很碎基本屬于偽脊線是噪聲引起的。我通常用屬性表里的長度字段做排序把幾百米以下的短線刪掉具體閾值看研究區(qū)尺度。這一步能顯著提高圖面質(zhì)量。5. 高頻問題排查這些坑我基本都踩過實(shí)際操作中總會碰到這樣那樣的問題。下面我把高頻問題整理成一個速查表方便你對照排查。問題現(xiàn)象可能原因解決辦法河網(wǎng)在山腳或平地處亂繞密密麻麻平行線平緩區(qū)域流量累積對微小地形起伏極其敏感先做坡度分析把坡度小于閾值的區(qū)域做掩膜排除或調(diào)大閾值山脊線偏離真實(shí)山脊常貼著山坡一側(cè)走反地形填洼過度限制Fill的Z limit檢查填洼前后差值只填淺洼地提取的山谷線不連續(xù)中間斷了很多段DEM含NoData、邊緣有無效值、洼地未填干凈先處理NoData檢查流向柵格是否有異常用“填洼Sink檢測”確認(rèn)線太細(xì)碎圖面全是毛刺分辨率太高、閾值太低、原始DEM噪聲大適當(dāng)提高閾值增加平滑容差刪除短小線對象山谷線和山脊線在鞍部交叉、串線鞍部地形復(fù)雜正反兩次提取都容易誤判將兩類線疊加到山體陰影上人工檢查必要時手動打斷編輯5.1 提取線斷裂八成問題出在DEM或流向“線不連續(xù)”是我在交流群里看到最多的一個問題。排查時我有一個固定套路先檢查原始DEM范圍和NoData邊界再看填洼前后的差值柵格。如果DEM邊緣某些位置完全沒有水流覆蓋提取結(jié)果就會在這些區(qū)域斷掉。還有一種情況是流向柵格中存在循環(huán)流向比如DEM里有局部平臺區(qū)域D8算法會把部分像元的水流方向指向自己導(dǎo)致流量累積出現(xiàn)斷點(diǎn)。遇到這種情況可以用柵格計算器把流向值為0或異常的像元單獨(dú)挑出來放大檢查平臺區(qū)域。如果DEM本身沒問題那問題大概率出在閾值上。閾值設(shè)置太高會把一些本來連續(xù)的細(xì)溝全部濾掉導(dǎo)致主干線斷斷續(xù)續(xù)。這時候不是去改DEM而是把閾值調(diào)低再配合線長篩選。5.2 閾值真的不是玄學(xué)是按需調(diào)試的很多人喜歡問“閾值到底填多少”我理解因?yàn)槟J(rèn)值往往效果不好。但這個問題真沒有標(biāo)準(zhǔn)答案。經(jīng)驗(yàn)參考是先用流量累積最大值的1%作為初值。如果研究區(qū)是典型的高山峽谷主溝與支溝級別差異大可以設(shè)置較小閾值比如0.2%到0.5%讓支溝也顯示出來如果是平緩丘陵閾值設(shè)置要大一些比如2%以上甚至更高。注意流量累積最大值不是固定不變的。填洼深度變了最大匯流量也會變所以換了DEM后一定要重新看統(tǒng)計值不能拿上一個項(xiàng)目的閾值直接套用。我的習(xí)慣是把閾值做成一個可調(diào)參數(shù)記錄在項(xiàng)目文檔里每次出圖時寫明用的什么閾值。這樣后續(xù)調(diào)整時非常高效而不是每次憑感覺亂填。如果你在做一個多期對比分析固定閾值還能讓不同期次的提取結(jié)果具備可比性。5.3 山谷線與山脊線交叉、串線的處理正向提取的山谷線和反向提取的山脊線在正常情況下應(yīng)該是相互交錯的山谷線“夾”在山脊線之間。但實(shí)際操作中常常會在鞍部、埡口這些地形復(fù)雜的位置出現(xiàn)交叉、串線。鞍部本身是山脊的最低點(diǎn)也是兩側(cè)山谷的源頭附近地形變化劇烈。正反兩次提取在鞍部都容易誤判正向提取可能把鞍部當(dāng)成溝谷的一部分反向提取又可能把鞍部誤認(rèn)為連續(xù)山脊。解決辦法沒有捷徑就是把山谷線和山脊線疊加到山體陰影或等高線底圖上人工檢查鞍部區(qū)域必要時手動打斷線、刪除錯誤線段。如果你需要做的是嚴(yán)格的拓?fù)浞治霰热绨焉郊咕€作為分水邊界來劃分子流域那還有另一個更嚴(yán)謹(jǐn)?shù)奶娲桨赣肁rcGIS水文分析里的Basin或Watershed工具直接生成分水嶺多邊形再把多邊形轉(zhuǎn)線得到分水線。這個結(jié)果比單純反向河網(wǎng)提取的山脊線更符合水文意義上的分水邊界。不過它的缺點(diǎn)是邊界會相對平滑細(xì)節(jié)不如反向河網(wǎng)豐富二者可以根據(jù)需求搭配使用。6. 最后補(bǔ)充幾條實(shí)際工作中的心得體會如果只用一句話總結(jié)這個流程我會說水文分析提取山谷山脊線本質(zhì)是兩次“逆向思維”的疊加第一次把山谷當(dāng)水溝第二次把山脊變成“假河谷”。掌握這個邏輯以后參數(shù)怎么調(diào)都不慌。再有幾條實(shí)操層面的體會DEM的預(yù)處理要占整個流程一半的功夫把NoData、投影、填洼處理好后面基本不會出大問題。輸出成果一定要疊加山體陰影做目視檢查這是最快、最直觀的質(zhì)量驗(yàn)證手段不要跳過。還有做項(xiàng)目報告時建議保留每一版閾值對應(yīng)的結(jié)果圖這樣方便論證參數(shù)選擇的合理性也方便返工調(diào)整。如果后續(xù)你需要在多個區(qū)域批量提取可以試試ArcGIS Model Builder把這套流程做成一個模型以后換DEM輸入直接跑。整個流程跑順之后我最大的感受是ArcGIS的水文分析工具其實(shí)是被低估的它不僅能做常規(guī)的河網(wǎng)提取稍微轉(zhuǎn)換一下思路就能解決很多地形分析里原本要手動勾繪的麻煩事。