基于DEM的大清河子流域劃分方法

張曉嬌 焦裕飛 劉佳 李旺林 李傳哲

摘?要：針對(duì)流域水循環(huán)模擬中的子流域劃分問(wèn)題，以大清河流域?yàn)槔紫然贏rcGIS軟件，以30 m分辨率的DEM作為數(shù)據(jù)源，進(jìn)行填洼處理、流向判定、匯流累計(jì)量計(jì)算等;然后采用河網(wǎng)密度法和水系分維法確定最佳集水面積閾值，提取出與實(shí)際最符合的河網(wǎng)，以此更準(zhǔn)確地提取流域水系;最后根據(jù)水流流向、匯流累計(jì)以及出水口的位置劃分子流域，并基于MATLAB程序?qū)ψ恿饔蚣捌鋬?nèi)部柵格進(jìn)行編碼。結(jié)果表明：當(dāng)集水面積閾值為270 km2時(shí)，提取的大清河流域河網(wǎng)水系與實(shí)際較為符合，且對(duì)子流域的劃分及編碼可以實(shí)現(xiàn)快速準(zhǔn)確定位每一個(gè)柵格，方便計(jì)算機(jī)的識(shí)別及演算。因此，以DEM為基礎(chǔ)實(shí)現(xiàn)大流域尺度下子流域劃分是可行的，對(duì)子流域的劃分及編碼可以為開(kāi)發(fā)構(gòu)建大尺度流域分布式水文模型提供技術(shù)支撐。

關(guān)鍵詞：DEM;水系提取;子流域劃分;大清河

中圖分類號(hào)：TV213.4文獻(xiàn)標(biāo)志碼：A

doi：10.3969/j.issn.1000-1379.2020.06.003

Study on Method of Sub-Basin Partition of Daqing River Based on DEM

ZHANG Xiaojiao1，2，JIAO Yufei2，LIU Jia2，LI Wanglin1，LI Chuanzhe2

（1. School of Water Conservancy and Environment，University of Jinan，Jinan 250022，China;

2. State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin，China

Institute of Water Resources and Hydropower Research，Beijing 100038，China）

Abstract：Aiming at the sub-basin partition in the basin water cycle simulation， taking Daqing River basin as an example， based on ArcGIS software and taking the DEM with the resolution of 30 m as data source for depression filling， determination of flow direction and calculation of cumulative amount of confluence， this paper adopted the river network density method and water system fractal method to determine the optimal catchment area threshold and extract the most practical river network， in order to extract the drainage system of Daqing River basin more accurately. Finally， the molecular watershed was divided according to the flow direction， confluence accumulation and the location of the outlet， and MATLAB program was used for the sub-watershed coding. The results show that when the threshold of catchment area is 270 km2， the extracted river network system of Daqing River basin is more consistent with the actual situation， and the partition and coding of sub-basins can achieve rapid and accurate positioning to each grid， which is convenient for computer identification and calculation. Therefore， it is feasible to realize subcatchment of large watershed scale based on DEM， and the coding of subcatchments can also provide technical support for the development and construction of distributed hydrological model of large watershed scale.

Key words： DEM;ditch extraction;subcatchments partition;Daqing River

子流域劃分是構(gòu)建分布式水文模型的基礎(chǔ)，分布式水文模型以子流域作為計(jì)算單元，可以清晰地了解單元之間的水文過(guò)程，并可很好地與傳統(tǒng)水文模型相結(jié)合，縮短模型開(kāi)發(fā)時(shí)間[1-2]。Jong-Gun K.等[3]研究認(rèn)為，子流域劃分對(duì)水文模擬結(jié)果的影響大于模型參數(shù)的影響，而且大尺度流域所劃分的子流域的數(shù)量不同使模擬精度存在不確定性。目前，基于數(shù)字高程模型（DEM）的流域地形特征提取方法已經(jīng)成為獲取流域信息的主要手段。OCallaghan等[4]基于水文學(xué)原理確定水流的路徑，并引入集水面積閾值的概念提取連續(xù)河網(wǎng)。傳統(tǒng)方法下作為子流域和河網(wǎng)劃分依據(jù)的集水面積閾值主要靠人為經(jīng)驗(yàn)設(shè)定，劃分結(jié)果與實(shí)際誤差較大。為了提高精度，學(xué)者們進(jìn)行了許多研究，目前確定最佳集水面積閾值的方法主要有河道平均坡降法[5]、河網(wǎng)密度法[6]、流域?qū)挾确植挤╗7]、水系分形法[8]等。 筆者針對(duì)大尺度流域的子流域劃分方法進(jìn)行研究，借助ArcGIS與MATLAB軟件，對(duì)大清河流域進(jìn)行河網(wǎng)提取、子流域劃分及編碼，以期為開(kāi)發(fā)構(gòu)建大清河流域分布式水文模型提供技術(shù)支撐。

1?子流域劃分方法

借助ArcGIS中的Hydrology工具箱以及MATLAB程序，對(duì)流域原始DEM進(jìn)行處理，步驟如下。

（1）地形預(yù)處理。原始的DEM中大多存在洼地或尖峰，會(huì)導(dǎo)致提取的水系出現(xiàn)斷裂或者水流方向錯(cuò)誤的現(xiàn)象[9]，因此需要對(duì)原始的DEM進(jìn)行填洼處理。

（2）提取水流流向。D8算法[4]分別用數(shù)字 1、2、4、8、16、32、64、128 來(lái)表示東、東南、南、西南、西、西北、北、東北8個(gè)方向（見(jiàn)圖1），利用該算法生成流域水流流向。

（3）計(jì)算匯流累計(jì)量。柵格匯流累計(jì)量的大小代表形成地表徑流的難易程度，數(shù)值越大，越易形成地表徑流。

（4）提取河網(wǎng)水系。通過(guò)設(shè)定一個(gè)合理的集水面積閾值，根據(jù)水流流向及匯流累計(jì)量提取河網(wǎng)水系[10]。采用河網(wǎng)密度法及水系分維法兩種方法來(lái)確定最佳集水面積閾值：①河網(wǎng)密度法，選取不同的集水面積閾值對(duì)流域水系進(jìn)行提取，計(jì)算不同閾值下的河網(wǎng)密度及河源密度，分別繪制集水面積閾值與河網(wǎng)密度、河源密度的關(guān)系曲線，曲線趨于平穩(wěn)的點(diǎn)即為最佳集水面積閾值;②水系分維法，采用盒維數(shù)法[11]，使用邊長(zhǎng)為r的正方形網(wǎng)格來(lái)覆蓋水系，確定覆蓋到水系的網(wǎng)格數(shù)N，分別取邊長(zhǎng)及網(wǎng)格數(shù)的對(duì)數(shù)進(jìn)行線性回歸，擬合直線的斜率即為水系的分維值，由此可以通過(guò)計(jì)算不同集水面積閾值下的水系分維值，與實(shí)際水系的分維值進(jìn)行比較，來(lái)確定最佳的集水面積閾值。

（5）河網(wǎng)分級(jí)。采用Strahler分級(jí)法對(duì)流域提取的水系進(jìn)行分級(jí)[12]，分級(jí)的原理是將河網(wǎng)中所有沒(méi)有支流匯入的河流作為第1級(jí)河流，兩個(gè)1級(jí)河流匯成第2級(jí)河流……依此類推，一直到河流出水口。在Strahler分級(jí)中，有且僅有同級(jí)別的兩條河流匯成一條河流時(shí)，河流的級(jí)別才會(huì)增大。

（6）確定子流域范圍。子流域的劃分是基于河流流向柵格和河流分級(jí)柵格，子流域內(nèi)匯流量最大的柵格為子流域的出水口，將所有能流向同一個(gè)子流域出水口的柵格劃在同一個(gè)子流域范圍內(nèi)[13]。

（7）子流域編碼。通過(guò)一定的編碼規(guī)則確定柵格之間的拓?fù)潢P(guān)系，以此來(lái)反映柵格的計(jì)算順序，便于計(jì)算機(jī)識(shí)別[14-15]?；跂鸥裰g的水流流向以及匯流關(guān)系，對(duì)每個(gè)柵格賦予一個(gè)5層編碼：第一層編碼A表示子流域的編碼;第二層編碼B表示該柵格在子流域中的匯流層數(shù)，子流域出水口柵格匯流層數(shù)為1，流入出水口柵格的柵格層數(shù)為2……依此類推;第三層編碼C表示匯流層內(nèi)的柵格順序，用來(lái)區(qū)分同一層內(nèi)的柵格;第四層編碼D表示流入柵格的匯流層;第五層編碼E用來(lái)區(qū)分河道與非河道柵格，其中河道柵格自出水口處從1開(kāi)始編碼，非河道柵格編碼為0。

2?應(yīng)用實(shí)例

2.1?研究區(qū)概況

大清河流域位于海河流域中部，流域面積4.3萬(wàn)km2，地跨京、津、冀、晉4?。ㄊ校?。大清河發(fā)源于太行山，分為白溝河水系和趙王河水系。該流域?qū)儆诘湫偷臏貛Ъ撅L(fēng)性氣候區(qū)，四季分明，降水量年際變化大、年內(nèi)分配不均，冬季和春季降雨少，夏季暴雨頻繁。

2.2?河網(wǎng)提取

借助ArcGIS軟件，對(duì)大清河流域原始DEM進(jìn)行填洼處理（見(jiàn)圖2），并利用D8算法提取水流流向（見(jiàn)圖3），計(jì)算匯流累計(jì)量。

河網(wǎng)提取的重點(diǎn)是通過(guò)設(shè)定一個(gè)最佳的集水面積閾值，使得提取的河網(wǎng)與實(shí)際河網(wǎng)水系相吻合。對(duì)于同一個(gè)流域，選取不同的集水面積閾值將得到不同的河網(wǎng)，閾值的選取也有很大的主觀性和隨意性，因此本文采用兩種方法確定最佳集水面積閾值。

（1）河網(wǎng)密度法。分別選取9、90、225、270、315、360、405、450 km2作為集水面積閾值，對(duì)大清河流域進(jìn)行水系提取，統(tǒng)計(jì)不同閾值下的河流長(zhǎng)度、河源數(shù)量，計(jì)算出相應(yīng)的河網(wǎng)密度與河源密度，結(jié)果見(jiàn)表1。

將集水面積閾值分別與河網(wǎng)密度、河源密度進(jìn)行趨勢(shì)線擬合（見(jiàn)圖4、圖5）。集水面積閾值與河網(wǎng)密度、河源密度都存在冪函數(shù)關(guān)系，擬合優(yōu)度R2都大于0.99。對(duì)所得到的冪函數(shù)求二階導(dǎo)數(shù)并再次進(jìn)行趨勢(shì)線擬合，將拐點(diǎn)作為最佳集水面積閾值（見(jiàn)圖6、圖7）。當(dāng)集水面積閾值為270 km2時(shí)，河網(wǎng)密度與河源密度二階導(dǎo)數(shù)都趨于0且之后不再變化，因此確定大清河流域水系提取的最佳集水面積閾值為270 km2。

（2）水系分形法。分別用邊長(zhǎng)r為100、500、1 000、2 000 m的網(wǎng)格對(duì)不同集水面積閾值下提取的水系進(jìn)行覆蓋，統(tǒng)計(jì)覆蓋水系的網(wǎng)格數(shù)目N，并用邊長(zhǎng)及網(wǎng)格數(shù)的對(duì)數(shù)之間的線性關(guān)系計(jì)算分維值，結(jié)果見(jiàn)表2。

分別用邊長(zhǎng)為100、500、1 000、2 000 m的網(wǎng)格在ArcGIS中對(duì)實(shí)際水系矢量圖進(jìn)行覆蓋，得到覆蓋水系的網(wǎng)格數(shù)，并用邊長(zhǎng)及網(wǎng)格數(shù)的對(duì)數(shù)之間的線性關(guān)系計(jì)算分維值，結(jié)果見(jiàn)表3。由表3可知，實(shí)際水系分維值為1.032，與集水面積閾值為270 km2時(shí)提取的水系分維值相同，因此確定最佳集水面積閾值為270 km2。

綜合上述兩種方法，可確定對(duì)大清河流域進(jìn)行水系提取的最佳集水面積閾值為270 km2。

2.3?河網(wǎng)分級(jí)

大清河流域河流分級(jí)見(jiàn)圖8。當(dāng)最佳集水面積閾值為270 km2時(shí)，大清河流域提取的河流一共有85條，其中：1級(jí)河流7條，占總河流數(shù)量的7.6%;2級(jí)河流10條，占總河流數(shù)量的10.7%;3級(jí)河流24條，占總河流數(shù)量的29.8%;4級(jí)河流44條，占總河流數(shù)量的51.9%。由此可見(jiàn)，大清河流域水系密布、發(fā)育較好。

2.4?子流域劃分

基于DEM提取的大清河流域面積為43 209 km2，與實(shí)際收集的數(shù)據(jù)基本吻合，可以認(rèn)為該DEM精度能夠滿足提取流域特征的要求。集水面積閾值為270 km2時(shí)，可將大清河流域劃分為56個(gè)子流域（見(jiàn)圖9），其中：子流域面積小于500 km2的占51%，在500～750 km2之間的占21%，在750～1 000 km2之間的占17%，大于1 000 km2的占11%。

2.5?子流域編碼

對(duì)子流域的編碼主要通過(guò)MATLAB程序來(lái)實(shí)現(xiàn)。圖10為大清河流域某一子流域內(nèi)所有柵格的編碼，以流域中某一柵格（18，3，5，3，150）為例，第一層編碼18表示第18個(gè)子流域，第二層編碼3為該柵格在第18個(gè)子流域中的第3層匯流層，第三層編碼5表示柵格為第3層匯流層中的第5個(gè)柵格，第四層編碼3表示流入該柵格的匯流所在層數(shù)為3，第五層編碼150表示第150條河道編碼。

3?結(jié)?論

本文采用河網(wǎng)密度法與水系分維法相互驗(yàn)證，確定大清河流域最佳集水面積閾值為270 km2，該條件下提取的水系與實(shí)際的水系吻合度很高。在此基礎(chǔ)上，應(yīng)用ArcGIS軟件對(duì)流域河網(wǎng)水系進(jìn)行分級(jí)并劃分子流域，然后借助MATLAB程序?qū)ψ恿饔蜻M(jìn)行編碼，可以實(shí)現(xiàn)對(duì)流域內(nèi)每一個(gè)柵格快速、精準(zhǔn)定位，有利于計(jì)算機(jī)的識(shí)別及演算。

通過(guò)對(duì)大清河流域水文特征進(jìn)行提取與研究，發(fā)現(xiàn)大尺度流域的子流域劃分關(guān)鍵在于確定最佳的集水面積閾值。由DEM數(shù)據(jù)提取的河網(wǎng)與實(shí)際誤差較小，可以滿足水文模擬的需要。通過(guò)基于水流流向與匯流計(jì)算的子流域劃分及編碼，可以充分了解分布式水文模型構(gòu)建過(guò)程中柵格之間水量交換等水文過(guò)程，為開(kāi)發(fā)適合大尺度流域的分布式水文模型奠定基礎(chǔ)。

參考文獻(xiàn)：

[1]?賈仰文，王浩，嚴(yán)登華.黑河流域水循環(huán)系統(tǒng)的分布式模擬（I）：模型開(kāi)發(fā)與驗(yàn)證[J].水利學(xué)報(bào)，2006，37（5）：534-542.

[2]?王中根，鄭紅星，劉昌明，等.黃河典型流域分布式水文模型及應(yīng)用研究[J].中國(guó)科學(xué)E輯：技術(shù)科學(xué)，2004，34（增刊1）：49-59.

[3]?JONG-GUN K，YOUNSHIK P，DONGSUN Y，et al. Development of a SWAT Patch for Better Estimation of Sediment Yield in Steep Sloping Watersheds[J]. Journal of the American Water Resources Association，2009，45（4）：963-972.

[4]?OCALLAGHAN J F. The Extraction of Drainage Networks from Digital Elevation Data[J].Computer Vision Graphics and Image Processing，1984，28（3）：323-344.

[5]?熊立華，郭生練.基于DEM的數(shù)字河網(wǎng)生成方法的探討[J].長(zhǎng)江科學(xué)院院報(bào)，2003，27（4）：14-17.

[6]?孔凡哲，李莉莉.利用DEM提取河網(wǎng)時(shí)集水面積閾值的確定[J].水電能源科學(xué)，2005，23（4）：65-67，93.

[7]?楊邦，任立良.集水面積閾值確定方法的比較研究[J].水電能源科學(xué)，2009，27（5）：11-14.

[8]?李麗. 分布式水文模型的匯流演算研究[D].南京：河海大學(xué)，2007：32-53.

[9]?徐新良，莊大方，賈紹鳳，等.GIS環(huán)境下基于DEM的中國(guó)流域自動(dòng)提取方法[J].長(zhǎng)江流域資源與環(huán)境，2004，13（4）：343-348.

[10]?宋向陽(yáng)，吳發(fā)啟，趙龍山，等.基于DEM的延河流域水文特征提取與分析[J].干旱地區(qū)農(nóng)業(yè)研究，2012，30（4）：200-206.

[11]?SKUBALSKA-RAFAJOWICZ E. A New Method of Estimation of the Box-Counting Dimension of Multivariate Objects Using Space-Filling Curves[J]. Nonlinear Analysis： Theory，Methods & Applications，2005，63（5-7）：e1281-e1287.

[12]?MELTON M A. A Derivation of Strahlers Channel-Ordering System[J]. Journal of Geology，1959，67（3）：345-346.

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

[14]?劉佳嘉，周祖昊，賈仰文，等.分布式水文模型子流域編碼方法對(duì)比分析[J].河海大學(xué)學(xué)報(bào)（自然科學(xué)版），2017，45（1）：22-29

[15]?王春，江嶺，陳泰生，等.基于Pfafstetter規(guī)則的流域編碼算法并行化方法[J].地球信息科學(xué)學(xué)報(bào)，2015，17（5）：556-561.

【責(zé)任編輯?張華興】