基于DEM的流域河網(wǎng)信息自動(dòng)提取算法

劉永和張萬(wàn)昌（ 河南理工大學(xué)資環(huán)學(xué)院，焦作 454000； 中國(guó)科學(xué)院遙感與數(shù)字地球研究所，北京 00094）

基于DEM的流域河網(wǎng)信息自動(dòng)提取算法

劉永和1張萬(wàn)昌2
（1 河南理工大學(xué)資環(huán)學(xué)院，焦作 454000；2 中國(guó)科學(xué)院遙感與數(shù)字地球研究所，北京 100094）

分布式陸面水文過(guò)程的模擬，除了經(jīng)典水文模型所必需的流域氣象、水文信息外，還需要研究流域的詳細(xì)地形、水系信息，方可實(shí)現(xiàn)流域內(nèi)產(chǎn)流和匯流的時(shí)空演算。以往通常需要借助商業(yè)軟件在研究區(qū)DEM上提取這些水系信息，不僅耗時(shí)、操作不便，還使得模型結(jié)構(gòu)松散。在自主研發(fā)的一種分布式水文模型的基礎(chǔ)上，開(kāi)發(fā)了一套用于提取水系信息的模塊，并與該水文模型以同一種語(yǔ)言緊密耦合，程序完全采用面向?qū)ο蟮姆绞阶灾骶帉?xiě)，大部分?jǐn)?shù)據(jù)交換均在內(nèi)存中直接進(jìn)行，而無(wú)需占用磁盤(pán)空間，運(yùn)行速度快，易于今后不斷完善和擴(kuò)展。詳細(xì)介紹了模塊中填洼與平坦區(qū)域處理、流向與累積流向矩陣生成、Strahler河道等級(jí)的確定、子流域生成、匯流次序和流程長(zhǎng)度信息的生成等多種流域河網(wǎng)信息提取的具體算法及程序?qū)崿F(xiàn)。本模型系統(tǒng)完全采用自主方式開(kāi)發(fā)，克服了以往使用商業(yè)軟件提取河網(wǎng)信息時(shí)的限制，使得分布式水文模式的流域模擬及分析功能更強(qiáng)大齊全和今后進(jìn)一步完善和擴(kuò)展。

分布式水文模型，數(shù)字高程模型，水系提取，耦合模塊

0　引言

模擬河流和土壤中水流的分布式水文模型被廣泛應(yīng)用于水利工程及環(huán)境影響評(píng)價(jià)[1-5]。分布式模型可用于估計(jì)與水流有關(guān)屬性的空間變化，如化學(xué)物質(zhì)遷移、土壤侵蝕、水分循環(huán)的陸面過(guò)程以及洪水狀況。由于流域水文狀況受天氣和氣候條件影響，人們期望在實(shí)現(xiàn)數(shù)值天氣預(yù)報(bào)和氣候變化預(yù)估的同時(shí)，能夠?qū)崿F(xiàn)流域的水文預(yù)報(bào)或未來(lái)水文條件預(yù)估。由于水流嚴(yán)格地受制于流域的地形特征，只有預(yù)先提取出流域的水系信息，才能實(shí)現(xiàn)分布式模擬計(jì)算地表徑流和壤中流的匯流。基于數(shù)字高程模型（DEM）提取流域水系是目前一種較為高效且可靠的手段。目前通過(guò)現(xiàn)代雷達(dá)測(cè)量手段已能獲得高精度的地面高程信息，如美國(guó)地質(zhì)調(diào)查局發(fā)布的SRTM DEM的精度已達(dá)到90m的水平分辨率，利用該數(shù)據(jù)能夠較為精準(zhǔn)地提取水系河網(wǎng)信息。因此，從DEM中由計(jì)算機(jī)自動(dòng)提取匯水網(wǎng)絡(luò)成為水文模擬研究中廣泛應(yīng)用的手段[6]。

通常由DEM提取的水系信息包括像元間的流向、河段以及對(duì)應(yīng)的子流域[7-8]。有關(guān)的研究已有很多文獻(xiàn)發(fā)表[3,7,9-13]?；诰匦我?guī)則格網(wǎng)的DEM是流域提取時(shí)應(yīng)用最多的數(shù)據(jù)[7]，但在一些文獻(xiàn)中也會(huì)基于不規(guī)則三角網(wǎng)（TIN）DEM[14-15]。本文主要指基于規(guī)則格網(wǎng)結(jié)構(gòu)的DEM。采用D8算法提取水系已經(jīng)接近成熟，具體細(xì)節(jié)見(jiàn)文獻(xiàn)[6, 16]。目前在一些成熟的商業(yè)軟件中已集成了有關(guān)基于D8方法的水系提取工具，如ArcGIS中的空間分析工具、ENVI軟件的插件River Tools等。然而基于商業(yè)軟件的水系提取工具其源代碼不公開(kāi)，人們無(wú)法知道其算法是否足夠高效且健壯，更重要的是商業(yè)軟件僅能作為一種數(shù)據(jù)處理工具，用戶(hù)無(wú)法對(duì)其處理過(guò)程和結(jié)果數(shù)據(jù)進(jìn)行自由操控，限制了它在自主開(kāi)發(fā)的分布式水文模型中的應(yīng)用。 分布式水文模型通常表現(xiàn)為一種執(zhí)行水文計(jì)算任務(wù)的計(jì)算機(jī)程序，由于缺少良好的人機(jī)交互界面而難以使用。目前一些重要的分布式水文模型（如SWAT）總是被集成于ArcGIS、GRASS等GIS軟件工具中，以充分借助GIS軟件所具有的DEM空間分析和可視化功能，增強(qiáng)其易用性。然而這種軟件集成方法對(duì)于分布式水文模型的研究者而言難度較大，影響了分布式水文模型的應(yīng)用。因此，開(kāi)發(fā)一套專(zhuān)門(mén)的水系信息提取工具，以解決當(dāng)分布式水文模型應(yīng)用中存在的問(wèn)題和不足具有一定的實(shí)際意義。

基于上述原因，我們針對(duì)分布式水文模型自主開(kāi)發(fā)了一套用于提供水系信息的耦合模塊。由于所建的分布式水文模型及水系提取模塊采用基于同一種語(yǔ)言（C#）的緊密耦合方式，大部分?jǐn)?shù)據(jù)交換均在內(nèi)存中直接進(jìn)行，而無(wú)需占用磁盤(pán)空間，且運(yùn)行速度快；程序完全采用面向?qū)ο蟮姆绞阶灾骶帉?xiě)，便于建立批量處理執(zhí)行任務(wù)，易于不斷完善和擴(kuò)展。

1　耦合系統(tǒng)的功能與耦合方案

1.1系統(tǒng)的功能

本文所述的流域河網(wǎng)信息自動(dòng)提取系統(tǒng)及其與分布式水文模型的耦合系統(tǒng)具有以下三類(lèi)功能：

1）類(lèi)似于GIS軟件中各種位圖和簡(jiǎn)單矢量圖顯示與漫游等可視化的功能。考慮到水文模型的主要任務(wù)是進(jìn)行流域水文模擬，本系統(tǒng)既具有可顯示單幅DEM或其他地理資料矩陣的功能，也可任意打開(kāi)并顯示用ARCGIS里所指定的ASCII二維數(shù)據(jù)格式或其他用ASCII碼矩陣方式存儲(chǔ)的二維數(shù)據(jù)。當(dāng)光標(biāo)在地圖上移動(dòng)時(shí)，系統(tǒng)可自動(dòng)計(jì)算并顯示當(dāng)前光標(biāo)處像元的值、行列號(hào)、經(jīng)緯度、UTM坐標(biāo)值。本系統(tǒng)可以作為查看任意兼容格式數(shù)據(jù)內(nèi)容的基本工具，如查看分布式水文模擬的降水、蒸發(fā)、徑流的空間分布狀況。為使顯示的DEM底圖或其他屬性圖像與矢量地圖相聯(lián)系，該系統(tǒng)還提供了設(shè)置和顯示點(diǎn)類(lèi)型和線(xiàn)類(lèi)型矢量圖的功能。矢量圖顯示的數(shù)據(jù)源為ARCGIS軟件導(dǎo)出的文本矢量格式。

2）自動(dòng)提取大量流域河網(wǎng)信息的功能。本系統(tǒng)中的流域河網(wǎng)提取模塊具有計(jì)算流向矩陣、累積流向矩陣、水系與河道等級(jí)矩陣、流域邊界、子流域劃分、匯流等級(jí)矩陣及匯流次序表、流程長(zhǎng)度矩陣和匯流時(shí)間矩陣等功能，為分布式水文模擬提供必要的信息。

3）分布式水文模型的模擬功能。本系統(tǒng)還具有按逐日和逐時(shí)進(jìn)行徑流和洪水模擬的功能，可按馬斯京根法實(shí)現(xiàn)河道洪水演算以及按等流時(shí)線(xiàn)滯時(shí)方法進(jìn)行匯流。

1.2系統(tǒng)的耦合方案

水文模型在運(yùn)行過(guò)程中所需要的氣象、土壤、植被、地形等空間分布數(shù)據(jù)均采用與所用DEM完全一致的分辨率。為方便使用該系統(tǒng)，在運(yùn)行過(guò)程中系統(tǒng)將以程序配置文件中所記錄的河流出口斷面所在的像元位置（行號(hào)和列號(hào)）作為運(yùn)行起始點(diǎn)，一次性完成所有河道信息提取過(guò)程。水文模擬過(guò)程中所需的所有河道信息以公用數(shù)據(jù)結(jié)構(gòu)變量的形式保存在內(nèi)存中，供水文模型中的產(chǎn)匯流模塊使用。

2　流域河網(wǎng)信息提取的算法

文中對(duì)流域河道信息的提取采用基于D8方法的思想。在每個(gè)DEM柵格上，待提取的信息包括流向、累積流向、水系與河道等級(jí)、流域邊界、子流域信息、匯流次序、匯流時(shí)間、流程長(zhǎng)度矩陣。在計(jì)算流向前，需要對(duì)DEM進(jìn)行填洼（depression filling）和處理平坦區(qū)域，以便能夠使DEM每個(gè)像元通過(guò)流向與主河道連接起來(lái)。

2.1填洼和平坦區(qū)域的處理

DEM中的洼地是指局部位置的像元高程低于外圍所有像元，使計(jì)算出的流向呈內(nèi)流現(xiàn)象，即流向鏈條不能與外圍的主河道相連接。DEM上的洼地可能是地表上真實(shí)的洼地，如小的湖泊，也可能是由于測(cè)量誤差以及分辨率較低所造成的假性洼地。在水文研究中，無(wú)論流域中真實(shí)存在的局部洼地，還是假洼地，都被視作流域中排水網(wǎng)絡(luò)的組成部分，即必須通過(guò)流向?qū)⑦@些區(qū)域與水系主河道相連接。填洼曾經(jīng)是DEM水系提取中最難解決的問(wèn)題，有關(guān)處理該問(wèn)題的算法較少，用傳統(tǒng)的一些算法（如Janson等[11]在1988年提出的方法）對(duì)一幅1000行×1000列的DEM進(jìn)行填洼需花費(fèi)數(shù)小時(shí)的時(shí)間，且因易出現(xiàn)錯(cuò)誤而使運(yùn)行不穩(wěn)定，但Planchon等2002年提出的快速算法已解決了該問(wèn)題[12]，本文使用了這種新算法，一般僅需數(shù)秒的時(shí)間即能完成，且運(yùn)行穩(wěn)定。

平坦區(qū)域中，由于所有像元的高程都相等，不存在任何流向，阻礙了水系的生成。目前常用做法是用不同程度的低于分辨率微小數(shù)值抬高平坦區(qū)域中像元的高程，使抬高的像元與其鄰接像元存在高程差，便于計(jì)算流向，并通過(guò)流向與外側(cè)水系相連接。

2.2水系生成與河道等級(jí)的確定

生成水系時(shí)需要完成流向計(jì)算、流向累積矩陣計(jì)算和確定河道等級(jí)。

在本系統(tǒng)中，我們采用一種新的鄰域位置編碼方法（圖1）。在這種編碼方法中，編碼值排列為三行。從左到右，第一行中的編碼設(shè)為0、1、2，第二行中的編碼設(shè)為3、4、5，第三行的編碼設(shè)為6、7、8。需要注意的是，在被設(shè)為4的3×3的編碼窗口中心的編碼不代表任何鄰域位置，而是代表了當(dāng)前指定像元的位置。由這些方向編碼很容易算出當(dāng)前指定像元的鄰域像元的行列號(hào)。假定當(dāng)前像元的位置為（i，j），該像元處的水流流向編碼為m，則其下游鄰接像元的行列號(hào)為（i+m/3-1，j+m%3-1），其中符號(hào)%表示取余運(yùn)算。如果一個(gè)流向編碼為n的鄰域像元的流向指向本像元，則其流向編碼應(yīng)為8-n?？梢?jiàn)，這種流向編碼方法非常有用，且很直接，這是因?yàn)楫?dāng)前感興趣像元的流出和流入像元的位置能夠高效計(jì)算出來(lái)。

本系統(tǒng)中采用單流向算法來(lái)計(jì)算流向，算法步驟是：（1）首先對(duì)DEM邊界上像素的流向設(shè)為指向DEM外部；（2）對(duì)除邊界外的所有內(nèi)部像素，比較其與周?chē)?個(gè)鄰接像元的高程差構(gòu)成的梯度，選擇高程低于當(dāng)前像素且梯度最大的鄰接像元所在的方向作為流向。這里的梯度是以當(dāng)前像元與鄰接像元的高程差除以像元距離D來(lái)計(jì)算。對(duì)上左下右四個(gè)緊鄰像元，像元距離D為1，而對(duì)角線(xiàn)上的四個(gè)像元，像元距離D為。

流向累積矩陣中每個(gè)像元上的累積流向值被定義為流過(guò)該像元的所有上游像元的總數(shù)，它代表了流過(guò)指定像元的上游匯水面積。常用的統(tǒng)計(jì)當(dāng)前像元累積流向的方法有兩種，一種是遞歸找出當(dāng)前像元的上游像元，統(tǒng)計(jì)這些找到的像元總數(shù)，即為累積流向值。另一種方法是逐個(gè)遍歷DEM中的每個(gè)像元，當(dāng)遍歷到某個(gè)像元時(shí)，該像元的所有下游像元的累積流向值加1。當(dāng)DEM中的像元全部遍歷一遍時(shí)，便得到了累積流向矩陣。

在得到流向矩陣和累積流向矩陣的基礎(chǔ)上，可進(jìn)一步對(duì)河流進(jìn)行基于Strahler分級(jí)方法[17]的水系分級(jí)，即把位于頂端的不再有分支的河流稱(chēng)為第1級(jí)河流，由兩個(gè)以上的同級(jí)河流匯合組成更高一級(jí)河流，如果多條河流同時(shí)匯聚，但最高級(jí)別的支流只有一條，則匯合后的河流級(jí)別與其最高級(jí)別的支流的級(jí)別相同。文獻(xiàn)[18]對(duì)Strahler河流分級(jí)方法進(jìn)行了較詳細(xì)的描述，采用從最上游的像素開(kāi)始向下游搜索，較為可行，但可能有一定的搜索難度。而本文則提出了另一種不需要從最上游像素開(kāi)始搜索的方法。具體步驟如下：

1）給定一個(gè)累積流向（匯流面積）閾值，認(rèn)為大于該閾值的像元位于河流上，不大于該閾值的不屬于河流；

2）在DEM中找出所有累積流向值（匯流面積）大于指定閾值的所有像元，即找出所有位于河網(wǎng)上的像元。為便于快速檢索，將這些河網(wǎng)像元放入一個(gè)以行號(hào)和列號(hào)為關(guān)鍵字的字典集合中；

3）對(duì)字典集合中的所有河網(wǎng)像元逐個(gè)遍歷，對(duì)找到的每個(gè)像元執(zhí)行第（4）～（8）步操作；

4）若當(dāng)前像元（用C1表示）已被分配子流域標(biāo)志，則返回第（3）步，即跳過(guò)當(dāng)前像元，對(duì)遍歷到的下一個(gè)像元進(jìn)行處理；

5）從當(dāng)前像元的8個(gè)鄰接像元中找出屬于C1上游的所有像元，若沒(méi)有找到這樣的像元，則將C1的河流級(jí)別設(shè)為1（屬初級(jí)河流）并返回第（3）步執(zhí)行下一趟循環(huán)。若找到了這樣的像元，執(zhí)行后續(xù)步驟；

6）對(duì)C1的上游像元，統(tǒng)計(jì)最大河流級(jí)別值m，以及擁有最大河流級(jí)別值的像元個(gè)數(shù)n，并判斷是否C1的全部上游像元已被標(biāo)記過(guò)河流級(jí)別，如果已被標(biāo)記則執(zhí)行后續(xù)步驟，否則返回第（3）步執(zhí)行下一趟循環(huán)；

7）如果在前面已找出的C1的上游像元個(gè)數(shù)為1時(shí)，則C1的河流級(jí)別設(shè)為與其上游像元相同的值；如果C1的上游像元多于1個(gè)時(shí)，且有兩個(gè)或更多像元級(jí)別具有相同的上游最大級(jí)別m時(shí)，C1的河流級(jí)別設(shè)為m+1；如果C1的上游像元多于1個(gè)，但只有一個(gè)上游像元具有最大級(jí)別m，則C1的河流級(jí)別設(shè)為m；

8）返回第（3）步；

9）如果第（2）～（8）步中沒(méi)有任何像元被標(biāo)記過(guò)，則退出；否則從字典集合中刪除已標(biāo)記過(guò)的像素記錄，并返回第（2）步繼續(xù)執(zhí)行；

10）結(jié)束后退出。

2.3子流域劃分

子流域一般是指流域面積超過(guò)一定面積的分支流域。文獻(xiàn)[19]中提出了一種子流域劃分方法，但算法不夠具體。文獻(xiàn)[20]中提出了具體的算法，但略顯復(fù)雜。本文采用思路與文獻(xiàn)[20]基本一樣，但采用下面的算法來(lái)完成：

1）生成一個(gè)空隊(duì)列Q，將代表河流出口處的像元設(shè)為子流域最低編號(hào)（如0或1），并將該像元的記錄加入隊(duì)列；

2）若隊(duì)列不空，從隊(duì)列中推出一個(gè)像元C1，執(zhí)行后續(xù)步驟，否則退出執(zhí)行；

3）從像元C1的8個(gè)鄰接方向中逐個(gè)查找C1的上游像元；

4）若找到C1的一個(gè)上游像元為C2，將C2加入隊(duì)列；

5）判斷C2的累積流向（流域面積）值大于給定閾值，且C2與C1之間的累積流向值之差也大于給定閾值的條件是否成立，若成立則將C2標(biāo)記為新子流域（其標(biāo)記值為C1的標(biāo)記值加1），若不成立，則C2的標(biāo)記與C1相同；

6）返回第（2）步繼續(xù)執(zhí)行。

上面算法第（5）步中的條件十分關(guān)鍵。C2作為C1的上游像元，其上游流域面積必須大于閾值才能成為一個(gè)子流域。另外，還需要保證C1存在除C2外的較大分支才能使C2成為新子流域的像元，這是因?yàn)?，如果除C2外的其他分支的匯流面積過(guò)小，小分支無(wú)法獨(dú)立成為子流域，而與C2、C1屬同一子流域。步驟（5）中的條件“C2與C1之間的累積流向值之差也大于給定閾值”正是為了保證C1存在除C2外的較大分支。

2.4匯流等級(jí)矩陣、匯流次序表、流程長(zhǎng)度與匯流時(shí)間矩陣的生成

在分布式水文模型的洪水模擬中，需要根據(jù)流域中每個(gè)像元間的流向關(guān)系從上游向下游采用馬斯京根方法進(jìn)行河道匯流。僅僅根據(jù)流向矩陣無(wú)法得出正確的匯流結(jié)果，仍需要根據(jù)匯流路徑的遠(yuǎn)近對(duì)所有柵格劃分匯流順序，從而得到整個(gè)流域或子流域的匯流等級(jí)矩陣和匯流次序表。對(duì)于同樣屬最低級(jí)別的河流像元，模擬時(shí)距離流域出口像元遠(yuǎn)的像元應(yīng)該先計(jì)算，距離流域出口近的像元應(yīng)后計(jì)算。匯流等級(jí)矩陣中的元素值為理論上根據(jù)到流域出口的距離來(lái)確定的匯流次序，等級(jí)值小的像元要先計(jì)算匯流，相同匯流等級(jí)的多個(gè)像元應(yīng)同時(shí)計(jì)算匯流。但從實(shí)際計(jì)算的角度來(lái)看，匯流等級(jí)相同的多個(gè)像元仍是分別來(lái)計(jì)算匯流的，但計(jì)算順序可是任意的。為方便，一般要導(dǎo)出匯流次序表，用以代替匯流等級(jí)矩陣。匯流次序表以匯流的實(shí)際計(jì)算次序存放流域中所有的像元的行列號(hào)，而不是采用矩陣形式存放。

流程長(zhǎng)度矩陣是與流域DEM完全對(duì)應(yīng)的二維矩陣（數(shù)組），其每個(gè)元素存放流域中對(duì)應(yīng)像元通過(guò)流向拓?fù)潢P(guān)系確定的通向流域出口的空間距離。匯流時(shí)間矩陣與流程長(zhǎng)度矩陣類(lèi)似，只是其每個(gè)元素存放的是對(duì)應(yīng)像元的“水流”流往流域出口所花費(fèi)的時(shí)間。匯流時(shí)間矩陣可以通過(guò)將流程長(zhǎng)度矩陣直接除以平均流速粗略得到，也可以通過(guò)將各個(gè)像元的水流流速與相鄰上下游像元間的梯度建立關(guān)系來(lái)精確導(dǎo)出。

匯流次序表可以通過(guò)匯流等級(jí)矩陣或流向矩陣導(dǎo)出。從流向矩陣導(dǎo)出的方法是從河流出口逆向上溯的方法搜索。由于流程長(zhǎng)度矩陣的計(jì)算以及通過(guò)梯度精確計(jì)算匯流時(shí)間矩陣時(shí)也需要采用同樣的逆向上溯方法來(lái)計(jì)算，可將三個(gè)計(jì)算任務(wù)用同一個(gè)算法來(lái)完成。算法仍借助隊(duì)列來(lái)實(shí)現(xiàn)，具體步驟如下：

1）將流域出口像元入隊(duì)列Q，并往匯流次序表中添加該像元，對(duì)應(yīng)該像元的流程長(zhǎng)度矩陣設(shè)為1（對(duì)匯流時(shí)間矩陣可采用類(lèi)似的方法賦值）；

2）當(dāng)隊(duì)列不空時(shí)，從中推出一個(gè)像元C，否則跳到第（6）步執(zhí)行；

3）從8個(gè)鄰接像元中查找屬于C的上游流入像元；

4）對(duì)找到的所有上游流入像元的流程長(zhǎng)度值進(jìn)行賦值，對(duì)斜對(duì)角方向的流入像元，流程長(zhǎng)度值為C的流程長(zhǎng)度加1.414，對(duì)上下左右四方向的流入像元，流程長(zhǎng)度值為C的流程長(zhǎng)度加1。用類(lèi)似的方法可對(duì)匯流時(shí)間矩陣進(jìn)行賦值。把找到的上游流入像元按任意順序添加入?yún)R流次序表的末尾；

5）返回第（2）步執(zhí)行；

6）將匯流次序表反轉(zhuǎn)；退出。

3　其他地面信息的生成

除與水系有關(guān)的信息外，分布式水文模型仍需要借助地面坡度和坡向信息來(lái)計(jì)算匯流速度和水分蒸散發(fā)量。計(jì)算坡度和坡向的方法可參見(jiàn)文獻(xiàn)[21]。

為使流域位置上DEM的地形狀況更加直觀，可生成漫反射光照明暗圖。在計(jì)算機(jī)圖形學(xué)中，漫反射光照采用Lambert余弦定理來(lái)實(shí)現(xiàn)[22]。Lambert余弦定理是根據(jù)入射光向量與反射面的法向來(lái)計(jì)算的，因此先要計(jì)算出DEM中的每個(gè)像元所處位置上的法向。由幾何學(xué)知道，空間中的平面法向?yàn)橥ㄟ^(guò)平面內(nèi)的任意兩個(gè)向量的向量積。顯然，只要已知平面中的三點(diǎn)，就能確定該平面的法向，因此三角形是最適于光照處理的面片。而規(guī)則格網(wǎng)DEM中每相鄰四個(gè)像元都能拼成兩個(gè)三角形（圖2）。需要指出的是，在拼三角形時(shí)需要保證三角形的頂點(diǎn)順序都為逆時(shí)針或都為順時(shí)針。設(shè)三角形的三個(gè)頂點(diǎn)分別為則三角形某兩邊上的矢量

當(dāng)每個(gè)DEM像元的三角形法向被算出后，再給定一束來(lái)自無(wú)窮遠(yuǎn)處的虛擬入射光，即可算出像元上的光強(qiáng)。

4　流域河網(wǎng)信息提取實(shí)例

本文截取了山東省境內(nèi)的沂河流域范圍的DEM數(shù)據(jù)，分辨率為240m，在此基礎(chǔ)上生成水文模擬時(shí)所用的各種與流域和河網(wǎng)有關(guān)的信息。圖3顯示的是系統(tǒng)的主窗體，其主要工作區(qū)中顯示了DEM圖像以及該區(qū)域內(nèi)的河流、縣界及水文測(cè)站的矢量信息。圖4顯示了DEM原圖及由其導(dǎo)出的Strahler河流等級(jí)、子流域分區(qū)、流程長(zhǎng)度、匯流次序等信息的圖像。圖4b顯示的由DEM導(dǎo)出的河網(wǎng)與現(xiàn)有的矢量主河道基本一致，這表明流域河網(wǎng)信息提取模塊中的填洼及計(jì)算累積流向的算法是較為可靠的。圖4c顯示了以沂河流域內(nèi)臨沂水文站為出口點(diǎn)的面積閾值超過(guò)2000像素的子流域劃分，由圖可見(jiàn)，劃分的子流域基本合理，但仍存在不同子流域面積差別較大的問(wèn)題。圖4d顯示了由DEM導(dǎo)出的沂河流域內(nèi)各像素至河流出口處的流程長(zhǎng)度，圖4e顯示了由DEM導(dǎo)出的各像元的匯流次序，這兩幅圖表明對(duì)應(yīng)的算法計(jì)算結(jié)果是正確的。圖4f顯示了由DEM直接生成的虛擬化光照明暗圖，該圖能夠反映研究區(qū)域地形起伏的狀況。

5　結(jié)語(yǔ)

本文介紹了在分布式水文模型中基于DEM提取流域河網(wǎng)信息的耦合模塊的算法設(shè)計(jì)以及耦合系統(tǒng)的功能和耦合方案，對(duì)水系生成和河道等級(jí)的確定、子流域生成、匯流次序及流程長(zhǎng)度信息的生成以及光照明暗圖的生成等算法進(jìn)行了詳細(xì)介紹。耦合模塊采用基于內(nèi)部數(shù)據(jù)結(jié)構(gòu)的緊密耦合方式與分布式水文模型耦合成一體化系統(tǒng)，對(duì)水文模型所需的各種流域河網(wǎng)信息采用自動(dòng)批量式生成，節(jié)省了對(duì)磁盤(pán)空間的占用，提高了運(yùn)行速度，同時(shí)也簡(jiǎn)化了對(duì)水文模擬時(shí)流域信息的準(zhǔn)備過(guò)程，能夠較大幅度的提高使用效率。耦合模塊中的各種算法結(jié)果基本合理或正確，運(yùn)行速度快，十分穩(wěn)定。

本水文模型系統(tǒng)的開(kāi)發(fā)完全采用自主方式，完全避免了使用商業(yè)GIS軟件在基于DEM準(zhǔn)備各種河網(wǎng)信息時(shí)的限制和不足，更有利于今后的進(jìn)一步完善和擴(kuò)展。

[1]Yang X J. Use of LIDAR elevation data to construct a highresolution digital terrain model for an estuarine marsh area. International Journal of Remote Sensing, 2005, 26(23): 5163-5166.

[2]Gumbo B, Munyamba N, Sithole G, et al. Coupling of digital elevation model and rainfall-runoff model in storm drainage network design. Physics and Chemistry of the Earth, 2002, 27(11-22): 755-764.

[3]Wang X H, Yin Z Y. A comparison of drainage networks derived from digital elevation models at two scales. Journal of Hydrology, 1998, 210(1-4): 221-241.

[4]Wolock D M, Mccabe G J. Differences in topographic characteristics computed from 100-and 1000-m resolution digital elevation model data. Hydrological Processes, 2000, 14(6): 987-1002.

[5]Dodov B A, Foufoula-Georgiou E. Floodplain morphometry extraction from a high-resolution digital elevation model: A simple algorithm for regional analysis studies. IEEE Geoscience and Remote Sensing Letters, 2006, 3(3): 410-413.

[6]Tribe A. Automated recognition of valley lines and drainage networks from grid digital elevation models - A review and a new method. Journal of Hydrology, 1992, 139(1-4): 263-293.

[7]Turcotte R, Fortin J P, Rousseau A N, et al. Determination of thedrainage structure of a watershed using a digital elevation model and a digital river and lake network.Journal of Hydrology, 2001, 240(3-4): 225-242.

[8]Paz A R, Collischonn W. River reach length and slope estimates for large-scale hydrological models based on a relatively hill highresolution digital elevation model. Journal of Hydrology, 2007, 343(3-4): 127-139.

[9]Colombo R, Vogt R V, Soille P, et al. Deriving river networks and catchments at the European scale from medium resolution digital elevation data. CATENA, 2007, 70(3): 296-305.

[10]Ahamed T, Rao K G, Murthy J. Automatic extraction of tank outlets in a sub-watershed using digital elevation models. Agricultural Water Management, 2002, 57(1): 1-10.

[11]Jenson S K, Domingue J O. Extracting topographic structure from digital elevation data for geographic information-system analysis. Photogrammetric Engineering and Remote Sensing, 1988, 54(11): 1593-1600.

[12]Planchon O, Darboux F. A fast, simple and versatile algorithm to fill the depressions of digital elevation models. CATENA, 2002, 46(2-3): 159-176.

[13]Jana R, Reshmidevi T V, Arun P S, et al. An enhanced technique in construction of the discrete drainage network from low-resolution spatial database. Computers & Geosciences, 2007, 33(6): 717-727.

[14]Tucker G E, Lancaster S T, Gasparini N M, et al. An objectoriented framework for distributed hydrologic and geomorphic modeling using triangulated irregular networks. Computers & Geosciences, 2001, 27(8): 959-973.

[15]劉學(xué)軍, 王永君, 任政, 等. 基于不規(guī)則三角網(wǎng)的河網(wǎng)提取算法.水利學(xué)報(bào), 2008, 39(1): 27-34.

[16]Mcmaster K J. Effects of digital elevation model resolution on derived stream network positions. Water Resources Research, 2002, 38(4): 13-1-13-8.

[17]Strahler A N. Quantitative analysis of watershed geomorphology. Transactions AGU, 1957, 38(6): 913-920.

[18]張渭軍, 王文科, 孔金玲, 等. 基于數(shù)字高程模型的水系自動(dòng)生成. 大地測(cè)量與地球動(dòng)力學(xué), 2006, 26(4): 41-44.

[19]沈中原, 李占斌, 李鵬, 等. 基于DEM的流域數(shù)字河網(wǎng)提取算法研究. 水資源與水工程學(xué)報(bào), 2009, 20(1): 20-28.

[20]葉愛(ài)中, 夏軍, 王綱勝, 等. 基于數(shù)字高程模型的河網(wǎng)提取及子流域生成. 水利學(xué)報(bào), 2005, 36(5): 531-537.

[21]湯國(guó)安, 劉學(xué)軍, 閭國(guó)年. 數(shù)字高程模型及地學(xué)分析的原理與方法. 北京: 科學(xué)出版社, 2005.

[22]唐澤圣, 周嘉玉, 李新友. 計(jì)算機(jī)圖形學(xué). 北京: 清華大學(xué)出版社, 1994.

A DEM Based Drainage Networks Extraction Module Seamlessly Integrated in a Distributed Hydrological Model

Liu Yonghe1, Zhang Wanchang2
(1 Institute of Resources and Environment, Henan Polytechnic University, Jiaozuo 454000 2 Institute of Remote Sensing and Digital Earth (RADI), Chinese Academy of Sciences, Beijing 100094)

For simulation of the hydrology in a watershed by distributed hydrological models, besides the meteorological and hydrological data needed, the detailed information such as topography and watershed networks is also necessary. Traditionally, the drainage networks were derived from digital elevation models (DEM) using commercial software, which is time consuming and inconvenient for operation in hydrological modeling. Based on a distributed hydrological model (DHM) developed by us, an automatic watershed information extraction module was developed which is able to be integrated into the DHM seamlessly using the same computer language of C#. Therefore, most of the data transferring can be finished in memory with no occupation on hard disks, so it is highly efficient when running. The algorithms for the watershed information extraction were introduced, such as the algorithms of removing of depression and flatten areas, generation of flow direction and accumulating flow direction, obtaining Strahler river order, dividing of sub watershed and calculation the sequence of flow order and water-drainage length. This system overcomes the limitations of using commercial software for extracting watershed information, and also favors the modeling system for convenience of modifying and extending the model easily in the future.

distributed hydrological models, digital elevation models, drainage networks extraction, coupled module

10.3969/j.issn.2095-1973.2015.02.009

2013年12月12日；

2014年3月11日

劉永和（1976—），Email: sucksis@163.com

資助信息：國(guó)家自然科學(xué)基金項(xiàng)目（41105074和41275108）；中科院數(shù)字地球重點(diǎn)實(shí)驗(yàn)室開(kāi)放基金項(xiàng)目（2011LDE010）；河南理工大學(xué)博士基金（B2011-038）