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

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

關鍵詞：DEM;水系提取;子流域劃分;大清河

中圖分類號：TV213.4文獻標志碼：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

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

1?子流域劃分方法

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

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

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

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

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

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

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

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

2?應用實例

2.1?研究區(qū)概況

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

2.2?河網(wǎng)提取

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

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

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

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

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

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

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

2.3?河網(wǎng)分級

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

2.4?子流域劃分

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

2.5?子流域編碼

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

3?結?論

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

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

參考文獻：

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

[2]?王中根，鄭紅星，劉昌明，等.黃河典型流域分布式水文模型及應用研究[J].中國科學E輯：技術科學，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].長江科學院院報，2003，27（4）：14-17.

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

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

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

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

[10]?宋向陽，吳發(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]?葉愛中，夏軍，王綱勝，等.基于數(shù)字高程模型的河網(wǎng)提取及子流域生成[J].水利學報，2005，36（5）：531-537.

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

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

【責任編輯?張華興】