基于數(shù)字高程模型的河網(wǎng)和流域邊界提取

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

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

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

1 研究區(qū)域與數(shù)據(jù)資料

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

覆蓋研究區(qū)域的DEM數(shù)據(jù)，用于水文特征的提取?？晒┻x擇的DEM數(shù)據(jù)源較多，如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]。可從網(wǎng)格分區(qū)中選取目標(biāo)區(qū)域的Geo TIFF文件。若目標(biāo)區(qū)域涵蓋幾個(gè)分區(qū)，則需做合并處理。本文選取美國(guó)奮進(jìn)號(hào)航天飛機(jī)的雷達(dá)地形測(cè)繪SRTM(Shuttle Radar Topography Mission)數(shù)據(jù)，目前最新為V4.1版本[8]，分辨率3″(約90 m)，高程垂直誤差小于16 m，水平分辨率高于現(xiàn)行全國(guó)智能網(wǎng)格預(yù)報(bào)(1～5 km)，提取數(shù)據(jù)精度能達(dá)到日常業(yè)務(wù)應(yīng)用要求。

2 主要算法原理及關(guān)鍵技術(shù)

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

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

圖2 基于DEM的數(shù)字化河網(wǎng)及流域邊界提取流程

2.1 DEM數(shù)據(jù)預(yù)處理

DEM數(shù)據(jù)預(yù)處理主要包括投影、拼接、裁剪與填洼。

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

一般用原始DEM數(shù)據(jù)計(jì)算流向(Flow Direction)，可以發(fā)現(xiàn)結(jié)果并不是期望的8個(gè)預(yù)定方位，而是被賦值1～255。因?yàn)橛袇R的存在，流向被賦予了所有可能流向的和。匯查找工具(Sink)可以計(jì)算得到所有匯。可以認(rèn)為所有匯，均不是自然存在[9]，而是數(shù)據(jù)原因?qū)е?，全部需要填充。填洼Fill，輸入原始DEM作為參數(shù)。可選參數(shù)z limit不填，即默認(rèn)全部填充。完成填洼運(yùn)算后，得到無(wú)凹陷DEM數(shù)據(jù)。

2.2 水流方向計(jì)算

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

圖3 方向編碼及流向

根據(jù)流向分析結(jié)果，可以計(jì)算每個(gè)單元格的上游匯流能力，如圖4所示。依據(jù)需要(如集水面積)設(shè)定匯流能力閾值，不低于該閾值的單元格標(biāo)記為河谷。不難發(fā)現(xiàn)，最大值處應(yīng)為該流域出水口。該方法簡(jiǎn)單，可直接產(chǎn)生連續(xù)的河網(wǎng)。

流向計(jì)算采用Flow Direction 工具。以填洼后的無(wú)凹陷DEM為輸入?yún)?shù)，可以選擇同時(shí)輸出drop raster，即D8算法中的距離權(quán)落差。這一步可以檢驗(yàn)填洼是否完全。正常情況，應(yīng)得到預(yù)定的8個(gè)方位編碼值，結(jié)果只能是1，2，4，……128，如若結(jié)果范圍為1～255，則證明填洼不完全，需要重復(fù)填洼步驟，直到結(jié)果為預(yù)定8方位編碼為止。

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

2.3 匯流累積量計(jì)算

采用Flow Accumulation工具，將流向結(jié)果作為參數(shù)輸入，權(quán)重參數(shù)默認(rèn)相同。因?yàn)楹泳W(wǎng)分級(jí)與流域劃分，一般按照一定的控制面積來(lái)確定，而像元大小均一。若通過(guò)徑流量的大小來(lái)確定，則需要設(shè)定每個(gè)網(wǎng)格的權(quán)重參數(shù)。因?yàn)槊總€(gè)網(wǎng)格的降水量和地質(zhì)地貌等條件都將影響該網(wǎng)格最后的徑流量，全區(qū)域并不均一。另外，此部分計(jì)算相對(duì)較慢，尤其當(dāng)DEM精度較高且研究區(qū)域較大時(shí)。

2.4 河網(wǎng)提取分級(jí)與矢量化

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

河網(wǎng)柵格矢量化，Stream to Feature可以將設(shè)定匯流累積量后的河網(wǎng)，編號(hào)后的河網(wǎng)連接以及分級(jí)河網(wǎng)等柵格數(shù)據(jù)轉(zhuǎn)化為矢量數(shù)據(jù)，以shapefile格式存儲(chǔ)，方便交換使用。

另外柵格矢量化Raster to Polyline工具，也可以矢量化河網(wǎng)，區(qū)別是前者主要針對(duì)河流網(wǎng)絡(luò)矢量化或者說(shuō)任何一類方向明確的代表線條網(wǎng)絡(luò)的柵格數(shù)據(jù)的矢量化，經(jīng)過(guò)優(yōu)化。表現(xiàn)在基于方向數(shù)據(jù)的輔助下，該算法將相同值的兩個(gè)相鄰線條特征矢量化為兩條平行線，而后者會(huì)將兩條線條重疊到一起，建議選擇Stream to Feature工具。

2.5 流域劃分與矢量化

流域提取工具有兩個(gè)，分別是流域盆地(Basin)和集水區(qū)(Watershed)劃分，主要區(qū)別是Basin從全局?jǐn)?shù)據(jù)范圍劃定流域，而Watershed可以通過(guò)設(shè)定出水口位置的方法，提取任意(子)流域。若無(wú)出水點(diǎn)位置柵格或者矢量數(shù)據(jù)，可以利用已經(jīng)生成的Stream Link作為匯水區(qū)的出水點(diǎn)，此方法因出水點(diǎn)較多，將得到當(dāng)前條件下最細(xì)小的子流域劃分。最后流域面可以通過(guò)矢量化工具Raster to Polygon轉(zhuǎn)為shapefile文件存儲(chǔ)。

2.6 提取檢驗(yàn)

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

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

3 分級(jí)河網(wǎng)及流域文件在氣象中的應(yīng)用

不同于單一圖形文件，利用數(shù)字化的河網(wǎng)和流域邊界信息，可制作GrADS掩膜mask文件。設(shè)定子流域編碼為1～n，在全域范圍內(nèi)按照氣象數(shù)值模式網(wǎng)格，將所有格點(diǎn)賦值其所在子流域編碼值，流域外則賦值為0?；诖颂囟╩ask文件，可方便利用GrADS函數(shù)(aave(maskout(YS,0.5-abs(mask.2-YsSubCode)),g))計(jì)算所有子流域面雨量，利用Draw String和Draw shp繪制填值、填色圖。圖5是歐洲中心細(xì)網(wǎng)格模式預(yù)報(bào)沅水流域未來(lái)12～84 h降水量的原始格點(diǎn)值與流域面雨量統(tǒng)計(jì)圖。

圖5 歐洲中心2017-05-19T08模式預(yù)報(bào)沅水流域未來(lái)12～84 h降水量(a 模式格點(diǎn)降水填色，b 子流域面雨量統(tǒng)計(jì))

可以發(fā)現(xiàn)預(yù)報(bào)降水呈北多南少的態(tài)勢(shì)，沅水流域大部分地區(qū)有中雨，南源龍頭江與南側(cè)支流渠水、巫水為小雨，各子流域統(tǒng)計(jì)結(jié)果與模式原始格點(diǎn)數(shù)據(jù)相符。子流域統(tǒng)計(jì)后的預(yù)報(bào)結(jié)論，匯流關(guān)系清晰，更利于防汛抗旱與水庫(kù)調(diào)度工作的開展。此方法可直接使用模式原始數(shù)據(jù)，充分利用了GrADS內(nèi)部函數(shù)，無(wú)需另行統(tǒng)計(jì)各子流域面雨量，且簡(jiǎn)便可行，符合氣象行業(yè)標(biāo)準(zhǔn)[14]。

4 結(jié)論與展望

本文以沅水流域?yàn)槔ㄟ^(guò)ArcGIS軟件，運(yùn)用D8算法，詳盡介紹了其中原理，展示了相關(guān)流程與實(shí)現(xiàn)步驟，并結(jié)合實(shí)踐中可能遇到的問題，重點(diǎn)提出了幾點(diǎn)細(xì)節(jié)。此方法，相對(duì)純手工繪制簡(jiǎn)便快捷。提取的河網(wǎng)及流域邊界文件，可直接應(yīng)用于MICAPS、SWAN等氣象業(yè)務(wù)系統(tǒng)。結(jié)合特定mask文件，還可利用氣象常用繪圖軟件GrADS的內(nèi)部函數(shù)，簡(jiǎn)便的應(yīng)用于流域面雨量預(yù)報(bào)的圖形展示。由于DEM精度等各方面原因，提取結(jié)果與實(shí)際情況存在一些細(xì)微出入，尤其在平原地區(qū)，因高程差小，相對(duì)誤差大，容易產(chǎn)生平行河流。充分利用高分辨率DEM數(shù)據(jù)結(jié)合衛(wèi)星影像的比對(duì)訂正，是提高水系流域提取精度的有效方法。