基于數字高程模型的河網和流域邊界提取

賈岸斌，佘高杰 ，盧 舟，蒯星龍

(1.常德市氣象局，湖南常德 415000；2.桃源縣氣象局，湖南桃源 415700)

河網和流域邊界的描述是水文氣象相關業務與研究的前提，是做好流域面雨量預報預測及預警服務等各項工作的前提條件[1]。傳統上(甚至現在)，這項工作一般是通過地形圖或等高線圖手工繪制完成的[2]。隨著遙感、GIS等行業及數字高程模型(Digital Elevation Modell，簡稱DEM)的發展，為數字化提取水文信息提供了可能。周成虎、杜云艷利用AVHRR影像，提出了基于水體光譜的水體自動提取識別模型[3]。DEM本身不包含河流、湖泊水庫、堤壩等信息，通過DEM直觀判讀流域水系及邊界較遙感影像困難，需要借助算法實現。目前最常用的是基于地形表面流水物理模擬分析算法中的D8算法。利用遙感影像識別水體和DEM提取流域邊界無論在精度上還是效率上都高于傳統地形圖手工勾繪方法。遙感影像具有獲取成本低、覆蓋率高、更新快、現勢性強等優點，但也存在云、積雪覆蓋、山體陰影等不利影響[2]。DEM在山地丘陵區和平均地形坡度不小于3°的區域所生成河網具有很高的可靠性，而在平均地形坡度小于3°的平坦區域河網生成中產生的虛擬河網與自然水系偏差較大[4]，但DEM還可以很方便地提取其它流域參數，仍然被廣為使用。本文以沅水流域為例，嘗試利用數字高程模型制作沅水主要河流的河網及流域邊界信息化文件，并將其應用于流域面雨量預報服務工作，希望能對氣象部門開展以流域為對象的相關預報服務及研究工作，比如針對河(湖)長的預報服務[5]等，起到一定的幫助。

1 研究區域與數據資料

沅水，又稱沅江，位于26°N～30°N、107°E～112°E之間，屬于長江流域洞庭湖水系，干流全長1 033 km，流域面積89 163 km2，跨湘黔渝鄂四省(市)。沅水自河源至黔城為上游，多深山幽谷，黔城至沅陵為中游，為丘陵地區，沅陵以下稱下游，桃源以下為沖積平原[6]，如圖1所示。

覆蓋研究區域的DEM數據，用于水文特征的提取。可供選擇的DEM數據源較多，如GTOPO30 (https://lta.cr.usgs.gov/GTOPO30)、SRTM DEM(http://srtm.csi.cgiar.org)、ASTER GDEM (http://gdem. ersdac.jspacesystems.or.jp/)、World DEM (http://www.astrium-geo.com) 等[7]。可從網格分區中選取目標區域的Geo TIFF文件。若目標區域涵蓋幾個分區，則需做合并處理。本文選取美國奮進號航天飛機的雷達地形測繪SRTM(Shuttle Radar Topography Mission)數據，目前最新為V4.1版本[8]，分辨率3″(約90 m)，高程垂直誤差小于16 m，水平分辨率高于現行全國智能網格預報(1～5 km)，提取數據精度能達到日常業務應用要求。

2 主要算法原理及關鍵技術

在GIS環境下，基于DEM的河網及流域邊界等水文信息的數字化提取流程如圖2所示，主要有DEM預處理，流向、匯流累積量計算等水文分析。最后將柵格數據矢量化，獲得GIS交換文件。

圖1 沅水流域地理位置及地形示意圖

圖2 基于DEM的數字化河網及流域邊界提取流程

2.1 DEM數據預處理

DEM數據預處理主要包括投影、拼接、裁剪與填洼。

一般DEM數據，只定義了地理坐標，無投影坐標，需要從Project工具，選擇WGS_1984_World_Mercator將原地理坐標GCS_WGS_1984投影到平面上。否則，雖不影響提取，但長度、面積屬性將無法使用。拼接(Mosaic To New Raster)之前，先導入一個DEM數據文件，查看相關屬性，方便拼接時填寫。主要有像素類型(16_BIT_SIGNED)及色帶數量(1)。DEM裁剪工具位于Spatial Analyst Tools-Extraction-Extract by Mask處。

一般用原始DEM數據計算流向(Flow Direction)，可以發現結果并不是期望的8個預定方位，而是被賦值1～255。因為有匯的存在，流向被賦予了所有可能流向的和。匯查找工具(Sink)可以計算得到所有匯。可以認為所有匯，均不是自然存在[9]，而是數據原因導致，全部需要填充。填洼Fill，輸入原始DEM作為參數。可選參數z limit不填，即默認全部填充。完成填洼運算后，得到無凹陷DEM數據。

2.2 水流方向計算

最先出現的流域特征提取方法是1975年Peuker等提出的識別谷點法[10]，但是該方法即便經過(Band和Donglas)改進仍然不能適用于大多數地貌類型[11]。目前應用最廣泛的是由奧克拉芳(O’Challaghan)和馬克(Mark)兩人于在1984年提出的無限坡面流累計方法—D8算法[12]。D8算法是最早出現的經典的單流向算法。該方法對自然狀態的水流方向進行了極大的概括，計算簡單，可操作性強。算法假設單個網格中的水流只有 8 種可能流向，即流入與之相鄰的8個網格中。它用最陡坡度法來確定水流的方向，在3×3的 DEM網格中，計算中心網格與各相鄰網格間的距離權落差(即網格中心點落差除以網格中心點之間的距離)，取距離權落差最大的網格為中心網格的流出網格，該方向即為中心網格的流向。流向分析如圖3所示。可以發現若定義網格距為1單位，則與其正交相鄰4個單元格中心距離也為1單位，與其對角相鄰4個單元格中心距離則為21/2，約1.414 21。

圖3 方向編碼及流向

根據流向分析結果，可以計算每個單元格的上游匯流能力，如圖4所示。依據需要(如集水面積)設定匯流能力閾值，不低于該閾值的單元格標記為河谷。不難發現，最大值處應為該流域出水口。該方法簡單，可直接產生連續的河網。

流向計算采用Flow Direction 工具。以填洼后的無凹陷DEM為輸入參數，可以選擇同時輸出drop raster，即D8算法中的距離權落差。這一步可以檢驗填洼是否完全。正常情況，應得到預定的8個方位編碼值，結果只能是1，2，4，……128，如若結果范圍為1～255，則證明填洼不完全，需要重復填洼步驟，直到結果為預定8方位編碼為止。

圖4 基于流向的匯流累積量分析

2.3 匯流累積量計算

采用Flow Accumulation工具，將流向結果作為參數輸入，權重參數默認相同。因為河網分級與流域劃分，一般按照一定的控制面積來確定，而像元大小均一。若通過徑流量的大小來確定，則需要設定每個網格的權重參數。因為每個網格的降水量和地質地貌等條件都將影響該網格最后的徑流量，全區域并不均一。另外，此部分計算相對較慢，尤其當DEM精度較高且研究區域較大時。

2.4 河網提取分級與矢量化

提取河網采用基于匯流累積量設定閾值的方法。閾值依據實際需要設定。例如網格分辨率為90 m×90 m，提取流域面積100 km2以上的河網數據，則匯流累計量需要不小于100×106/(90 m×90 m)個網格。通過Raster Calculator工具，SetNull("FlowAcc"<12345,1)，將河谷柵格設置為1，非河谷柵格設置為空(也可使用Con 工具)，得到河網柵格。河網連接Stream Link和河網分級Stream Order分別為每段河流分配唯一編號與標定相應級別。

河網柵格矢量化，Stream to Feature可以將設定匯流累積量后的河網，編號后的河網連接以及分級河網等柵格數據轉化為矢量數據，以shapefile格式存儲，方便交換使用。

另外柵格矢量化Raster to Polyline工具，也可以矢量化河網，區別是前者主要針對河流網絡矢量化或者說任何一類方向明確的代表線條網絡的柵格數據的矢量化，經過優化。表現在基于方向數據的輔助下，該算法將相同值的兩個相鄰線條特征矢量化為兩條平行線，而后者會將兩條線條重疊到一起，建議選擇Stream to Feature工具。

2.5 流域劃分與矢量化

流域提取工具有兩個，分別是流域盆地(Basin)和集水區(Watershed)劃分，主要區別是Basin從全局數據范圍劃定流域，而Watershed可以通過設定出水口位置的方法，提取任意(子)流域。若無出水點位置柵格或者矢量數據，可以利用已經生成的Stream Link作為匯水區的出水點，此方法因出水點較多，將得到當前條件下最細小的子流域劃分。最后流域面可以通過矢量化工具Raster to Polygon轉為shapefile文件存儲。

2.6 提取檢驗

為了檢驗提取效果，將沅水流域面積[13]與其對應提取面積對比，如表1，可以發現沅水流域面積誤差率在0.71%，其主要支流提取誤差除武水、酉水偏差4%左右略大外，其他均在±1%內。

表1 沅水流域面積提取誤差

3 分級河網及流域文件在氣象中的應用

不同于單一圖形文件，利用數字化的河網和流域邊界信息，可制作GrADS掩膜mask文件。設定子流域編碼為1～n，在全域范圍內按照氣象數值模式網格，將所有格點賦值其所在子流域編碼值，流域外則賦值為0。基于此特定mask文件，可方便利用GrADS函數(aave(maskout(YS,0.5-abs(mask.2-YsSubCode)),g))計算所有子流域面雨量，利用Draw String和Draw shp繪制填值、填色圖。圖5是歐洲中心細網格模式預報沅水流域未來12～84 h降水量的原始格點值與流域面雨量統計圖。

圖5 歐洲中心2017-05-19T08模式預報沅水流域未來12～84 h降水量(a 模式格點降水填色，b 子流域面雨量統計)

可以發現預報降水呈北多南少的態勢，沅水流域大部分地區有中雨，南源龍頭江與南側支流渠水、巫水為小雨，各子流域統計結果與模式原始格點數據相符。子流域統計后的預報結論，匯流關系清晰，更利于防汛抗旱與水庫調度工作的開展。此方法可直接使用模式原始數據，充分利用了GrADS內部函數，無需另行統計各子流域面雨量，且簡便可行，符合氣象行業標準[14]。

4 結論與展望

本文以沅水流域為例，通過ArcGIS軟件，運用D8算法，詳盡介紹了其中原理，展示了相關流程與實現步驟，并結合實踐中可能遇到的問題，重點提出了幾點細節。此方法，相對純手工繪制簡便快捷。提取的河網及流域邊界文件，可直接應用于MICAPS、SWAN等氣象業務系統。結合特定mask文件，還可利用氣象常用繪圖軟件GrADS的內部函數，簡便的應用于流域面雨量預報的圖形展示。由于DEM精度等各方面原因，提取結果與實際情況存在一些細微出入，尤其在平原地區，因高程差小，相對誤差大，容易產生平行河流。充分利用高分辨率DEM數據結合衛星影像的比對訂正，是提高水系流域提取精度的有效方法。