基于DEM的丹江口水源區(qū)治理區(qū)小流域劃分研究

郭文慧 ，于 泳 ，李 璐 ，楊 偉 ，袁修猛

(1.湖北省水利水電科學(xué)研究院 湖北省水土保持工程技術(shù)研究中心，湖北 武漢 430070；2.湖北省水利廳，湖北 武漢 430070)

隨著計(jì)算機(jī)技術(shù)和地理信息系統(tǒng)的發(fā)展，數(shù)字高程模型DEM(Digital Elevation Model)中蘊(yùn)涵著大量的地形、地貌信息，可以有效地反映區(qū)域的基本地形空間分布規(guī)律與地貌特征，通過分析提取河網(wǎng)、流域界線甚至提取流域內(nèi)的地形屬性都是切實(shí)可行的。隨著數(shù)字水文的興起和分布式水文模型研究的發(fā)展，利用DEM獲取水系和子流域特征，進(jìn)而進(jìn)行水文模擬的方法越來越廣泛地為水文學(xué)者所采用。

本研究以1∶5萬DEM數(shù)字高程模型為基礎(chǔ)，在ArcMap軟件的水文模塊(Hydrology)環(huán)境下，探索了流域水系自動(dòng)提取及流域劃分的方法，并依據(jù)《小流域劃分及編碼規(guī)范》(SL 653—2013)的相關(guān)要求對(duì)子流域進(jìn)行了人工修正、歸并，最后完成了符合規(guī)范要求的適宜面積的小流域劃分。

1 研究區(qū)概況

以湖北省境內(nèi)國(guó)家水土保持重點(diǎn)治理區(qū)劃中的丹江口水源區(qū)治理區(qū)(以下簡(jiǎn)稱“治理區(qū)”)為研究對(duì)象，包括十堰市的張灣區(qū)、茅箭區(qū)、鄖西縣、鄖縣、丹江口市、竹山縣、竹溪縣、房縣共8個(gè)縣(市、區(qū))。治理區(qū)位于秦嶺東西向構(gòu)造體系的南部邊緣，地形的主要特點(diǎn)是高差大、坡度陡、切割深，總的地勢(shì)是西北高、東南低、北陡南緩，漢江沿線地貌峽谷和盆地相間。在我國(guó)氣候區(qū)劃上位于南北氣候過渡地帶的秦巴山區(qū)，北亞熱帶季風(fēng)氣候，多年平均氣溫為15.9 ℃，年均降水量約900 mm，年均水面蒸發(fā)量860 mm左右。土壤以山地黃棕壤和黃褐土為主，山地棕壤土次之，并有紫色土發(fā)育。

丹江口水庫是我國(guó)水資源調(diào)整的重大工程——南水北調(diào)中線工程水源地，1973年建成，設(shè)計(jì)蓄水位157 m，水庫面積745 km2，庫容174.5億m3(其中漢江庫區(qū)占53.9%、丹江庫區(qū)占46.1%)，漢江、丹江回水長(zhǎng)度分別為174和83 km。南水北調(diào)中線工程實(shí)施后，丹江口大壩加高至176.6 m，設(shè)計(jì)蓄水位170 m，總庫容290.5億m3，水庫面積1 050 km2，回水長(zhǎng)度漢江194 km、丹江93 km。漢江庫區(qū)的主要入庫支流有金錢河、將軍河、天河、曲遠(yuǎn)河、堵河、神定河、泗河、芝河、浪河等,丹江庫區(qū)的主要入庫支流有丹江、淇河、老灌河等。

研究區(qū)DEM情況見圖1。

圖1 研究區(qū)DEM情況

2 數(shù)據(jù)來源

(1)DEM數(shù)據(jù)。研究區(qū)1∶5萬DEM數(shù)據(jù)(柵格大小為25 m×25 m，見圖1)，投影坐標(biāo)系為WGS84坐標(biāo)系，UTM投影，6度分帶的49號(hào)帶。

(2)水系數(shù)據(jù)。 研究區(qū)5 km以上的河流水系圖。

(3)其他數(shù)據(jù)。研究區(qū)1∶5萬各縣級(jí)行政區(qū)劃矢量數(shù)據(jù)。

3 研究方法

由DEM提取流域的數(shù)字特征，包括DEM預(yù)處理、流向分析、 匯流累積分析及河網(wǎng)提取、流域邊界生成等過程。

(1)DEM預(yù)處理。從 DEM中自動(dòng)提取流域的前提是DEM中不能存在洼地和尖峰，否則在計(jì)算水流方向時(shí)會(huì)出現(xiàn)逆流的現(xiàn)象，給水流線和流域界線的確定帶來困難，為此首先要將所有地形改造成由斜坡構(gòu)成[1]。

(2)流向分析。在填洼后的DEM中，每個(gè)柵格單元都有一個(gè)可以定義的水流方向值。水流方向的確定有單流向法和多流向法兩種。單流向法包括D8算法、Rh08算法、Lea算法和D∞算法等，多流向法則包括MFD算法、DAEMOM算法和Ding算法等[2]。本研究采用目前較常用的D8算法,該算法是假定每個(gè)單元格的水流方向可能流向8個(gè)相鄰的單元格，水流沿最陡坡度的方向流動(dòng)，計(jì)算出中心單元格與周圍8個(gè)單元格的坡度即可確定中心單元格的水流方向[3]。

(3)匯流累積分析。從每個(gè)柵格單元出發(fā)依次掃描流向矩陣，沿水流方向追蹤到DEM邊界，當(dāng)整個(gè)水流方向矩陣掃描完畢，就可以得到流域匯流能力的柵格分布圖，匯流柵格上每個(gè)單元的值代表上游匯流區(qū)內(nèi)流入該單元的上游柵格單元格的總數(shù)，值較大者可視為河谷，值等于零則是較高的地方，可能為流域的分水嶺。根據(jù)匯流累積柵格分布圖可以方便地提取流域的各種特征參數(shù)[4]。

(4)河網(wǎng)提取。首先設(shè)定最小支流上游集水區(qū)的面積閾值，以上游集水區(qū)面積大于閾值面積的網(wǎng)格點(diǎn)作為該支流的起始點(diǎn)，流域內(nèi)集水面積超過該閾值的網(wǎng)格點(diǎn)定義為水道[5]。將匯流累積柵格上所有大于或等于最小集水面積閾值的柵格提取出來，即得到了柵格河網(wǎng)，然后將其轉(zhuǎn)化為矢量河網(wǎng)。

(5)流域邊界生成。流域又稱集水區(qū)域，是指流經(jīng)其中的水流和其他物質(zhì)從一個(gè)公共的出水口排出，從而形成一個(gè)集中的排水區(qū)域。集水流域的生成需要先確定該集水區(qū)域的最低點(diǎn)，然后結(jié)合水流方向數(shù)據(jù)，分析搜索出該出水點(diǎn)上游所有流過該出水點(diǎn)的柵格，直到所有的該集水區(qū)的柵格都確定了位置，也就是搜索到流域的邊界或分水嶺的位置。

(6)小流域歸并。本研究以《小流域劃分及編碼規(guī)范》(SL 653—2013)為依據(jù)進(jìn)行小流域歸并工作。小流域面積原則上控制在30～50 km2，特殊情況不宜小于3 km2或大于100 km2。

4 流域特征提取

4.1 無洼地DEM的生成

4.1.1 洼深計(jì)算

首先借助Hydrology工具集中的Flow Direction工具獲取水流方向數(shù)據(jù)。ArcGIS中的水流方向是利用D8算法，即最大距離權(quán)落差(最大坡降法)來確定的。

運(yùn)用Sink工具計(jì)算獲得洼地區(qū)域，并用Watershed工具生成洼地的貢獻(xiàn)區(qū)域，分別用Zonal Statistic和Zonal Fill工具計(jì)算得到洼地所形成的貢獻(xiàn)區(qū)域最低高程和洼地出水口高程。洼地深度即為出水口高程與貢獻(xiàn)區(qū)域最低高程之差。通過計(jì)算，研究區(qū)DEM最大洼地深度為88 m。

4.1.2 填洼閾值確定

分別將30、50、70、90 m作為填充深度值，在Hydrology工具中集中利用Fill工具將其處理，并對(duì)填洼后的DEM進(jìn)行匯流累積分析。對(duì)比分析匯流累積圖水系特征信息與已配準(zhǔn)的河流水系圖，確定最佳填洼閾值。

比較發(fā)現(xiàn)，30 m填洼結(jié)果坡面洼地有所減少，50 m填洼結(jié)果較為理想，70 m填洼結(jié)果與50 m填洼結(jié)果基本相同，90 m填洼結(jié)果所生成的匯流累積圖中的水系特征與河流矢量圖差異較大。為避免過度填洼，選擇50 m為填洼閾值。

4.2 河網(wǎng)的自動(dòng)提取

為使獲取的河網(wǎng)更加準(zhǔn)確，即水系的分級(jí)和水系的長(zhǎng)度更加符合實(shí)際，在閾值設(shè)定時(shí)進(jìn)行了一系列試驗(yàn)，分別選取5 km2(8 000個(gè)柵格)、10 km2(16 000個(gè)柵格)、20 km2(32 000個(gè)柵格)作為最小匯流面積閾值，河網(wǎng)提取結(jié)果依次見圖2、3、4。

將不同閾值條件下的河網(wǎng)提取結(jié)果及所提取河網(wǎng)與水系圖比較發(fā)現(xiàn)：5 km2匯流面積條件下提取的河網(wǎng)偏密，且河流長(zhǎng)度與實(shí)際水系一致或稍短；10 km2匯流面積條件下提取的河網(wǎng)稍密，且河流長(zhǎng)度較實(shí)際水系稍短；20 km2匯流面積下提取的河網(wǎng)稀疏，且河流長(zhǎng)度較實(shí)際水系短。

圖2 5 km2匯流面積的河網(wǎng)提取結(jié)果

圖3 10 km2匯流面積的河網(wǎng)提取結(jié)果

圖4 20 km2匯流面積的河網(wǎng)提取結(jié)果

4.3 流域邊界生成

對(duì)比分析表明，不同閾值條件下所生成的小流域數(shù)量與河道數(shù)一致，且一一對(duì)應(yīng)。經(jīng)統(tǒng)計(jì)，最小匯流閾值為5 km2時(shí)生成子流域2 923個(gè)，最大子流域面積為63.83 km2，50 km2以上子流域9個(gè)；最小匯流閾值為10 km2時(shí)生成子流域1 261個(gè)，最大子流域面積為131.04 km2，50 km2以上子流域49個(gè)；最小匯流閾值為20 km2時(shí)生成子流域615個(gè)，最大子流域面積為194.86 km2，50 km2以上子流域161個(gè)。

5 小流域邊界處理與劃分結(jié)果

5.1 子流域邊界處理與歸并

5.1.1 子流域邊界整合

根據(jù)《小流域劃分及編碼規(guī)范》(SL 653—2013)中的相關(guān)規(guī)定：小流域面積原則上控制在30～50 km2，特殊情況下不宜小于3 km2或大于100 km2。為提高數(shù)據(jù)處理效率并保證小流域邊界的準(zhǔn)確性，本研究選擇子流域面積不大于50 km2，且子流域數(shù)量相對(duì)較少的劃分結(jié)果用于小流域的編碼及命名。為此，選擇最小匯流閾值為10 km2時(shí)所提取生成的子流域邊界數(shù)據(jù)，對(duì)于大于50 km2的子流域，用最小匯流閾值為5 km2的提取結(jié)果進(jìn)行數(shù)據(jù)替換，并將小于0.1 km2的子流域與周邊子流域自動(dòng)合并，最終得到替換合并更新后的數(shù)據(jù)(子流域1 742個(gè)，大于50 km2的子流域2個(gè))用于小流域的修正、分割及歸并。

5.1.2 子流域邊界修正

由于基于DEM的河網(wǎng)提取是采用最大坡降的方法，因此在平地區(qū)域(例如谷底等)水流方向是隨機(jī)的，就很容易生成平行狀的河流等錯(cuò)誤形狀，出現(xiàn)大量復(fù)雜鑲嵌的流域格局(如圖5)， 此時(shí)需要利用Editor工具欄進(jìn)行手工編輯修正。

圖5 復(fù)雜鑲嵌的流域邊界示例

5.1.3 子流域邊界分割

首先在屬性表中添加“亞單元標(biāo)識(shí)碼”，選出跨縣的小流域，賦值為1，其他未跨縣小流域賦值為0；然后用湖北省分縣行政邊界矢量數(shù)據(jù)對(duì)跨縣子流域進(jìn)行分割。

5.1.4 子流域邊界歸并

依照《小流域劃分及編碼規(guī)范》(SL 653—2013)中的原則對(duì)子流域邊界進(jìn)行歸并，受小流域面積的限制，歸并具有一定的主觀性?，F(xiàn)將這次歸并研究過程中所遇到的特殊問題處理如下：

(1)被縣界分割形成面積較小的子流域。對(duì)此類子流域研究其與周邊子流域的關(guān)系，符合歸并原則的與周邊子流域歸并得到面積適合的小流域；對(duì)不符合歸并原則，且面積小于3 km2的子流域，為方便流域管理，考慮與周邊小流域合并；面積大于3 km2的子流域可單獨(dú)作為小流域。

(2)丹江口水庫及周邊區(qū)域。丹江口水庫周邊區(qū)域及庫底地形復(fù)雜，基于DEM自動(dòng)提取的子流域邊界也較為復(fù)雜，較多子流域邊界跨越水庫水面，與對(duì)岸一起形成區(qū)間型小流域，不利于小流域的管理?；诖耍紤]將170 m淹沒線以內(nèi)區(qū)域(包括丹江口水庫水面及庫周170 m淹沒線以內(nèi)區(qū)域)劃分為丹江口水庫型流域，將周邊的子流域邊界依據(jù)水系匯流關(guān)系及行政邊界進(jìn)行歸并。

(3)盲區(qū)。在省界邊緣，會(huì)出現(xiàn)少數(shù)面積較小的流入其他省份的水系，他們與本省內(nèi)周邊相鄰子流域的匯流關(guān)系較遠(yuǎn)。這部分子流域多小于5 km2，分析認(rèn)為造成這種結(jié)果一般有兩種原因：①實(shí)際地形地貌復(fù)雜；②省級(jí)邊界與DEM影像之間在地理位置上存在誤差。由于提取出來的水系與湖北省河流分布圖一致性較好，因此可排除第二種原因?？紤]到便于流域管理，依據(jù)周邊水系的匯流關(guān)系進(jìn)行歸并，形成類似于坡面小流域的類型(如圖6)。

圖6 特殊情況示例

5.2 小流域劃分結(jié)果

本研究在丹江口水源區(qū)治理區(qū)約2.30萬km2的區(qū)域內(nèi)，提取小流域624個(gè)(丹江口水庫型流域除外，該流域面積為389.04 km2)，最小的小流域面積為3.41 km2，最大的小流域面積為59.42 km2。其中：<30 km2的小流域有122個(gè)，面積約為0.25萬km2，占該區(qū)總面積的10.87%；30～50 km2的小流域有455個(gè)，面積約為1.80萬km2，占該區(qū)總面積的78.26%；≥50 km2的小流域有47個(gè)，面積約為0.25萬km2，占該區(qū)總面積的10.87%。研究區(qū)小流域劃分結(jié)果見圖7。

圖7 研究區(qū)小流域劃分結(jié)果

6 結(jié)論與討論

利用1∶5萬比例尺DEM數(shù)據(jù)提取的河網(wǎng)與實(shí)際河網(wǎng)的比較結(jié)果一致性較好；在提取水系圖時(shí)，由于水系為連續(xù)的匯流最低點(diǎn)相連接而形成的線狀圖，因此水系只能反映河流的走向，而不能反映河流的寬度。自動(dòng)提取的子流域的邊界與分水嶺基本吻合，流域特征明顯，但由于河底或庫底地形的復(fù)雜性，出現(xiàn)了大量復(fù)雜鑲嵌的流域格局，存在較多子流域需要人工修正，而人工修正又存在主觀性。

小流域是水土流失治理工作的基本單元。本研究對(duì)典型區(qū)域的小流域劃分技術(shù)進(jìn)行了研究，為湖北省開展類似區(qū)域小流域劃分工作提供了方法。

[參考文獻(xiàn)]

[1] 陳加兵，勵(lì)惠國(guó)，鄭達(dá)賢，等．基于DEM的福建省小流域劃分研究[J]．地球信息科學(xué)，2007，9(2):74-77.

[2] GONG J Y，XIE J B.Extraction of drainage networks from large terrain datasets using high through put computing[J]. Comp-uters & Geosciences，2009，35:337-346.

[3] 曹玲玲，張秋文．基于SRTM的數(shù)字河網(wǎng)提取及其應(yīng)用[J]．人民長(zhǎng)江，2007，38(8):150-152.

[4] 王俊啓．基于1∶1萬地形圖的數(shù)字流域劃分研究——以浙江省安吉縣為例[D].杭州：浙江大學(xué)，2013：1-55.

[5] 秦福來，王曉燕，王麗華，等．基于DEM的流域特征提取及應(yīng)用[J].河北師范大學(xué)學(xué)報(bào)(自然科學(xué)版)，2006,30(1):109-112.