劉 芝

(四川省資陽(yáng)市雁江區(qū)水務(wù)局，四川 資陽(yáng)，641300)

1 概述

在水庫(kù)洪水預(yù)報(bào)中，傳統(tǒng)的預(yù)報(bào)方法是繪制水庫(kù)洪水預(yù)報(bào)綜合查算圖[1]，但該查算圖需要根據(jù)本流域大量的實(shí)測(cè)降雨流量資料進(jìn)行配線，可移用性較差，必須做到一庫(kù)一圖，而對(duì)于實(shí)測(cè)資料缺乏的小型水庫(kù)，只能借用鄰近流域特性參數(shù)或地區(qū)綜合經(jīng)驗(yàn)參數(shù)繪制查算圖，預(yù)報(bào)結(jié)果的可信度不高。隨著社會(huì)經(jīng)濟(jì)的發(fā)展和社會(huì)管理水平的提高，尤其在“智慧地球”、“智慧城市”、“智慧水利”的發(fā)展框架下，國(guó)家在逐漸筑牢大江大河防洪體系的同時(shí)，對(duì)水庫(kù)洪水預(yù)報(bào)能力的要求也在不斷提高。研究和揭示水庫(kù)洪水與流域地形地貌的關(guān)系，建立揭示流域水文響應(yīng)與地形地貌關(guān)系的地貌瞬時(shí)單位線(GIUH)及具有物理基礎(chǔ)的分布式水文模型，是破解小型水庫(kù)洪水預(yù)報(bào)難題的必然出路。隨著計(jì)算機(jī)技術(shù)和地理信息技術(shù)的發(fā)展，尤其是數(shù)字高程模型(DEM)的廣泛應(yīng)用，提取流域河網(wǎng)及地形地貌信息成為可便捷實(shí)現(xiàn)的過(guò)程。地貌瞬時(shí)單位線(GIUH)和分布式水文模型[2]在小型水庫(kù)洪水預(yù)報(bào)中的推廣和應(yīng)用，成為未來(lái)的必然趨勢(shì)。

2 R-V GIUH基本理論

Rodriguez-Iturbe和Valdes于1979年首先提出了地貌瞬時(shí)單位線(GIUH)的概念，并將瞬時(shí)均勻降落在流域上的凈雨看作由無(wú)數(shù)的水質(zhì)點(diǎn)組成，用概率論的方法來(lái)研究流域上水質(zhì)點(diǎn)的運(yùn)動(dòng)，當(dāng)組成凈雨的水質(zhì)點(diǎn)間呈弱相互作用時(shí)，流域瞬時(shí)單位線等價(jià)于水質(zhì)點(diǎn)持留時(shí)間的概率密度函數(shù)，這就是經(jīng)典的R-V GIUH理論。

2.1 河流分級(jí)

R-V GIUH理論是基于Horton-Strahler河流分級(jí)法提出的。Horton-Strahler河流分級(jí)的具體方法為：將所有河網(wǎng)弧段中沒有支流的河網(wǎng)弧段定為第1級(jí)；兩個(gè)1級(jí)河網(wǎng)弧段匯成的河網(wǎng)弧段定為第2級(jí)；兩個(gè)2級(jí)河網(wǎng)弧段匯成的河網(wǎng)弧段定為第3級(jí)；如此下去，當(dāng)且僅當(dāng)同級(jí)別的兩條河網(wǎng)弧段匯成的河網(wǎng)弧段級(jí)別增加一級(jí)，直到河網(wǎng)出水口。

2.2 Horton地貌參數(shù)

R-V GIUH理論將瞬時(shí)單位線參數(shù)同流域地貌參數(shù)結(jié)合，提供了流域響應(yīng)函數(shù)的地貌學(xué)解釋。自然界的水系在發(fā)育過(guò)程中遵循一定的規(guī)律。由于受地質(zhì)構(gòu)造和自然環(huán)境的影響，河網(wǎng)形態(tài)在平面上表現(xiàn)為有規(guī)律的排列并在不同尺度上具有一定的自相似性。Horton通過(guò)大量研究發(fā)現(xiàn)[3]，任何一個(gè)流域內(nèi)，各級(jí)河道的數(shù)目、平均長(zhǎng)度、平均集水面積、平均縱比降與河道的Horton-Strahler河流分級(jí)級(jí)別之間存在幾何級(jí)數(shù)關(guān)系，并據(jù)此提出了Horton地貌定律，即河數(shù)定律、河長(zhǎng)定律、面積定律、比降定律。而Horton地貌定律中的河數(shù)率RB、河長(zhǎng)率RL、面積率RA是推求R-V地貌瞬時(shí)單位線時(shí)重要的地貌參數(shù)。

Horton地貌參數(shù)(河數(shù)率RB、河長(zhǎng)率RL、面積率RA)的計(jì)算公式如下：

RB=Ni-1/Nii=2，3，…，Ω

一般來(lái)說(shuō)，不同級(jí)別的河流的河數(shù)率、河長(zhǎng)率、面積率并不相同，應(yīng)用時(shí)可取各級(jí)河流的平均值。且大量研究成果顯示：天然河網(wǎng)的Horton地貌參數(shù)存在一個(gè)較窄的變化范圍，河數(shù)率范圍為3～5，河長(zhǎng)率范圍為1.5～3.5，面積率范圍為3～6。

2.3 R-V GIUH基本公式推導(dǎo)

基于Horton-Strahler河流分級(jí)方法，把水質(zhì)點(diǎn)到達(dá)的某一級(jí)河流或坡面稱作狀態(tài)xi，把水質(zhì)點(diǎn)從注入點(diǎn)到出口斷面順序經(jīng)過(guò)的各狀態(tài)稱作路徑si，則水質(zhì)點(diǎn)沿路徑si流至出口斷面的匯流時(shí)間為Tsi=Tx1+Tx2+…+TxΩ。根據(jù)流域瞬時(shí)單位線等價(jià)于水質(zhì)點(diǎn)持留時(shí)間的概率密度函數(shù)，可導(dǎo)出流域地貌瞬時(shí)單位線基本表達(dá)式[4]：

μ(t)=fB(t)=∑siSfx1(t)×fx2(t)×…×fxΩ(t)·p(si)

(1)

式中：fx1(t)，fx2(t)，…，fxΩ(t)表示水質(zhì)點(diǎn)在當(dāng)前狀態(tài)下的等待時(shí)間概率密度函數(shù)；p(si)表示水質(zhì)點(diǎn)取路徑si的概率；S表示所有路徑si的集合，有S={s1，s2，……}；Ω表示河流最高級(jí)別；符號(hào)×表示卷積相乘。

式(1)中的路徑概率p(si)有：

p(si)=θx1px1x2px2x3…pxΩ-1xΩ

(2)

式中：θxi為水質(zhì)點(diǎn)處于初始狀態(tài)的概率，簡(jiǎn)稱初始概率；pxixj(i

(3)

(4)

式(1)中等待時(shí)間概率密度函數(shù)fxi(t)的確定，實(shí)質(zhì)上就是確定水質(zhì)點(diǎn)流速的概率密度函數(shù)。常用下列單參數(shù)指數(shù)型函數(shù)來(lái)進(jìn)行擬合：

(5)

式中：Ki為水質(zhì)點(diǎn)在狀態(tài)xi的平均持留時(shí)間，等于狀態(tài)的空間尺度與水質(zhì)點(diǎn)在該狀態(tài)的平均流速的比值，可用地形地貌學(xué)方法確定。

將R-V GIUH基本表達(dá)式導(dǎo)出為數(shù)學(xué)計(jì)算公式的過(guò)程比較復(fù)雜，且當(dāng)河流最高級(jí)別Ω不同時(shí)，導(dǎo)出的計(jì)算公式形式亦不相同，使用時(shí)需要分河網(wǎng)級(jí)別予以導(dǎo)出。R-V GIUH理論原作者給出了三級(jí)河網(wǎng)的R-V GIUH計(jì)算公式，由文康[6]等于1988年給出了四級(jí)河網(wǎng)和五級(jí)河網(wǎng)的R-V GIUH計(jì)算公式。

2.4 流域平均流速的推求

(6)

式中：v為流域平均流速；L為主河道長(zhǎng)度；τ為流域匯流時(shí)間，可用下列公式推求：

(7)

(8)

式中：ψ為洪峰徑流系數(shù)；F為流域集水面積；L為主河道長(zhǎng)度；J為河道平均比降；S為暴雨雨力；n為暴雨公式指數(shù)；m為匯流參數(shù)。以上參數(shù)均能利用《四川省水文手冊(cè)》進(jìn)行查算。

3 基于DEM提取流域河網(wǎng)和Horton地貌參數(shù)

基于DEM提取流域河網(wǎng)和地貌參數(shù)可在ArcGIS軟件中實(shí)現(xiàn)。ArcGIS軟件是美國(guó)環(huán)境系統(tǒng)研究所ESRI開發(fā)的一款無(wú)縫擴(kuò)展GIS軟件，具有強(qiáng)大的地理信息分析處理能力。ArcGIS軟件中的Hydrology工具集能夠方便地提取河網(wǎng)結(jié)構(gòu)等大量水文信息。

3.1 DEM的獲取

數(shù)字高程模型(DEM)是表示地面高程隨空間變化的數(shù)學(xué)模型，是對(duì)地形地貌的一種離散的數(shù)字表達(dá)。DEM的獲取方式主要有三種：一是通過(guò)互聯(lián)網(wǎng)上既有的DEM數(shù)據(jù)庫(kù)直接下載，常用的DEM產(chǎn)品有ASTERGDEM、SRTM、GLS2005、TANDEM等，分辨率多為90m和30m；二是通過(guò)現(xiàn)有的地形等高線圖生成DEM，即將現(xiàn)有地形等高線圖數(shù)字矢量化后，進(jìn)行TIN建模，進(jìn)而生成DEM；三是根據(jù)實(shí)測(cè)數(shù)據(jù)源生成DEM，DEM數(shù)據(jù)源的獲取方式包括地面測(cè)量、數(shù)字?jǐn)z影、激光雷達(dá)、立體遙感、GPS等。普通的研究和應(yīng)用采用實(shí)測(cè)數(shù)據(jù)源生成DEM的可能性不高，通常采用第一、二種DEM獲取方式，且以互聯(lián)網(wǎng)下載方式最為便捷。

3.2 DEM流域位置定位

互聯(lián)網(wǎng)下載獲得的DEM是分區(qū)塊的，在區(qū)塊DEM數(shù)據(jù)上定位到具體研究區(qū)域的位置，需要借助于行政邊界矢量數(shù)據(jù)。在ArcGIS軟件中打開DEM數(shù)據(jù)，并將下載到的行政邊界矢量數(shù)據(jù)加載到DEM上(需注意坐標(biāo)系統(tǒng)的一致性)，然后使用ArcToolbox中的Extract by Mask工具進(jìn)行剪切，即得到行政邊界區(qū)域內(nèi)的DEM數(shù)據(jù)。將對(duì)應(yīng)行政區(qū)域的水系圖通過(guò)空間校正疊加到DEM數(shù)據(jù)上，可定位出具體流域的位置，再沿集水區(qū)域邊界剪切，便可得到具體流域范圍的DEM數(shù)據(jù)。

3.3 DEM洼地填充處理

受DEM分辨率的限制，以及DEM生產(chǎn)過(guò)程中的系統(tǒng)誤差，原始DEM表面存在著一些凹陷或尖峰的區(qū)域，這些是DEM數(shù)據(jù)的缺陷，會(huì)導(dǎo)致不合理的甚至是錯(cuò)誤的水流方向計(jì)算結(jié)果。但并非所有的洼地都是數(shù)據(jù)誤差造成的，也有可能是真實(shí)存在的地形。故需要進(jìn)行洼地計(jì)算，計(jì)算出洼地的位置、區(qū)域、深度，以此判斷哪些洼地是數(shù)據(jù)誤差，哪些洼地是真實(shí)地形。具體計(jì)算步驟如下：①用Hydrology工具集下的FlowDirection工具計(jì)算出水流方向；②依據(jù)水流方向，用Sink工具計(jì)算出洼地位置；③用Watershed工具計(jì)算出洼地貢獻(xiàn)區(qū)域；④用Zonal工具集下的ZonalStatistic工具計(jì)算出洼地貢獻(xiàn)區(qū)域最低高程；⑤用Zonal Fill工具計(jì)算出洼地貢獻(xiàn)區(qū)域出口最低高程；⑥用MapAlgebra工具集下的Raster Calculator工具計(jì)算出洼地深度，洼地深度=洼地貢獻(xiàn)區(qū)域出口最低高程-洼地貢獻(xiàn)區(qū)域最低高程。依據(jù)計(jì)算出的洼地深度設(shè)置合理的洼地填充閾值，用Hydrology工具集下的Fill工具對(duì)洼地進(jìn)行填平處理，得到接近于真實(shí)地形的無(wú)洼地DEM。

3.4 河網(wǎng)提取

河網(wǎng)是水在地形地貌上的產(chǎn)物，同時(shí)水的匯集與流動(dòng)又受到河網(wǎng)的約束。河網(wǎng)特征反映一個(gè)流域的綜合水文特征。常用的河網(wǎng)提取方法為連續(xù)河網(wǎng)生成法，即采用D8法按最陡坡度確定DEM中每個(gè)網(wǎng)格單元的水流方向，根據(jù)水流方向計(jì)算出每個(gè)網(wǎng)格單元的上游集水面積(即匯流累積量)，設(shè)置一個(gè)集水面積閾值，以該閾值作為河道認(rèn)定的門檻，由匯流累積量大于該閾值的網(wǎng)格單元組成河網(wǎng)。河網(wǎng)提取具體步驟如下：①用Hydrology工具集下的Flow Direction工具計(jì)算出無(wú)洼地DEM的水流方向；②依據(jù)無(wú)洼地DEM水流方向，用Flow Accumulation工具計(jì)算出匯流累積量；③設(shè)置適當(dāng)?shù)募娣e閾值，用Map Algebra工具集下的RasterCalculator工具生成柵格河網(wǎng)；④用Hydrology工具集下的Stream to Feature工具將柵格河網(wǎng)轉(zhuǎn)換為矢量河網(wǎng)。

3.5 流域劃分

流域又稱集水區(qū)域，是指由分水線所包圍的河流集水區(qū)，流經(jīng)其中的水流從一個(gè)公共的出水口流出。流域的邊界可依據(jù)DEM的數(shù)字地形分析確定。流域劃分具體步驟如下：①依據(jù)提取出的流域河網(wǎng)，用Hydrology工具集下的Stream Link工具計(jì)算出河網(wǎng)結(jié)點(diǎn)；②依據(jù)無(wú)洼地DEM水流方向，用Basin工具劃分出流域盆地；③依據(jù)河網(wǎng)結(jié)點(diǎn)計(jì)算結(jié)果，用Snap Pour Point工具計(jì)算出匯水區(qū)出水口；④用Watershed工具生成集水流域；⑤用Conversion Tools工具集下的Raster to Polygon工具對(duì)柵格集水流域進(jìn)行矢量化，得到矢量集水流域。

3.6 Horton地貌參數(shù)的提取

Horton地貌定律是以Horton-Strahler河流分級(jí)為基礎(chǔ)的。河網(wǎng)的分級(jí)包含重要的水文信息，可用以研究水流的運(yùn)行和匯流模式等。河網(wǎng)分級(jí)具體步驟為：①根據(jù)提取的柵格河網(wǎng)，用Hydrology工具集下的Stream Order工具，選擇Strahler分級(jí)法，完成河網(wǎng)的分級(jí)；②用Stream to Feature工具對(duì)柵格分級(jí)河網(wǎng)進(jìn)行矢量化，得到矢量分級(jí)河網(wǎng)。

Horton地貌參數(shù)的提取步驟具體如下：①打開矢量分級(jí)河網(wǎng)的屬性列表，添加列，用Calculate Geometry工具計(jì)算出河段長(zhǎng)度，導(dǎo)出含河段長(zhǎng)度的屬性列表，便可獲得每一級(jí)別河段的條數(shù)和長(zhǎng)度等信息；②打開矢量集水流域的屬性列表，添加列，用Calculate Geometry工具計(jì)算出集水流域面積，導(dǎo)出含集水流域面積的屬性列表，便可獲得每一級(jí)河流的集水面積等信息；③依據(jù)Horton地貌參數(shù)計(jì)算公式，計(jì)算得到河數(shù)率RB、河長(zhǎng)率RL、面積率RA等地貌參數(shù)值。

4 應(yīng)用實(shí)例

4.1 研究區(qū)域概況

雁江區(qū)位于四川盆地中部，全區(qū)共有小型水庫(kù)78座，水庫(kù)流域面積介于0.1km2～450km2之間，普遍沒有降雨和流量觀測(cè)站，受預(yù)報(bào)方法和管理水平的限制，其水庫(kù)洪水預(yù)報(bào)工作基本上是空白的。本次研究以高板橋水庫(kù)流域作為研究對(duì)象，進(jìn)行流域河網(wǎng)和Horton地貌參數(shù)的提取，并建立水庫(kù)流域的R-V地貌瞬時(shí)單位線。高板橋水庫(kù)位于雁江區(qū)東峰鎮(zhèn)境內(nèi)，屬沱江水系，地理坐標(biāo)為E104°53′18.3″、N30°3′50.3″，總庫(kù)容402萬(wàn)m3，壩址以上集雨面積33.77km2，主河道長(zhǎng)13.355km，河道平均坡比降1.317‰，是一座以防洪為主、兼顧農(nóng)業(yè)灌溉等綜合效益的小(1)型水庫(kù)。

4.2 流域河網(wǎng)和地貌參數(shù)的提取

研究使用的DEM為從地理空間數(shù)據(jù)云免費(fèi)下載的ASTER GDEM V2 30m分辨率DEM，雁江區(qū)的行政邊界矢量數(shù)據(jù)為從全國(guó)地理信息資源目錄服務(wù)系統(tǒng)免費(fèi)下載獲得，流域定位所用的雁江區(qū)水系圖比例尺為1∶85000。通過(guò)圖層疊加和剪切，可分別獲得雁江區(qū)和高板橋水庫(kù)流域的DEM，見圖1。

(a)雁江區(qū)DEM (b)高板橋水庫(kù)流域DEM

按照洼地計(jì)算方法對(duì)原始DEM進(jìn)行洼地計(jì)算，得到洼地的位置、區(qū)域和深度。經(jīng)計(jì)算，高板橋水庫(kù)流域共有洼地21處，洼地深度最大值25m，最小值1m。由于本研究使用的原始DEM分辨率只有30m，故結(jié)合實(shí)際地形考慮，認(rèn)為以上洼地均是由于數(shù)據(jù)誤差產(chǎn)生的偽洼地，故將所有洼地予以填平處理，得到研究區(qū)域的無(wú)洼地DEM。

設(shè)置適當(dāng)?shù)募娣e閾值，利用ArcGIS軟件中的Hydrology工具集進(jìn)行流域河網(wǎng)提取，得到高板橋水庫(kù)流域的流域河網(wǎng)，見圖2。

(a)雁江區(qū)流域河網(wǎng) (b)高板橋水庫(kù)流域河網(wǎng)

對(duì)提取出的高板橋水庫(kù)流域河網(wǎng)進(jìn)行Strahler分級(jí)運(yùn)算，獲得具有河流分級(jí)屬性的流域河網(wǎng)數(shù)據(jù)，并以此計(jì)算出各級(jí)河流的條數(shù)、長(zhǎng)度、集水面積。根據(jù)Horton地貌參數(shù)計(jì)算公式，可計(jì)算得到高板橋水庫(kù)流域的河數(shù)率、河長(zhǎng)率、面積率等地貌參數(shù)，見表1。

表1 高板橋水庫(kù)流域Horton地貌參數(shù)

可見，基于DEM提取出的高板橋水庫(kù)流域河網(wǎng)共分成4級(jí)，總集水面積42.0783km2，與從萬(wàn)分之一地形圖上量算到的水庫(kù)壩址以上集雨面積33.77km2存在一定差異。高板橋水庫(kù)流域的河數(shù)率為4.3667、河長(zhǎng)率為2.3821、面積率為4.9665，均在自然河流Horton地貌參數(shù)的取值范圍內(nèi)，可用于進(jìn)行流域R-V地貌瞬時(shí)單位線的推求。

4.3 R-V地貌瞬時(shí)單位線推求

根據(jù)《四川省水文手冊(cè)》用推理公式求得高板橋水庫(kù)流域的平均匯流時(shí)間約為8h，進(jìn)而求得高板橋水庫(kù)流域平均流速約1.67m/s。

將提取的高板橋水庫(kù)流域地貌參數(shù)和計(jì)算的流域平均流速帶入四級(jí)河網(wǎng)流域地貌瞬時(shí)單位線公式進(jìn)行計(jì)算，發(fā)現(xiàn)初始概率θx4<0，這是不合理的，此處采用文康[6]等推薦的處理方式，即當(dāng)不滿足RB/RA<0.8這個(gè)約束條件，初始概率出現(xiàn)負(fù)值時(shí)，直接以不包括低一級(jí)河道集水面積在內(nèi)的各級(jí)河網(wǎng)集水面積與流域總面積之比表示初始概率。

最后，可得出高板橋水庫(kù)流域的地貌瞬時(shí)單位線為：

μ(t)=-0.1671e-9.2492t-20.5336e-4.5271t+6.5354e-0.9442t+17.2233e-3.058t

將上述地貌瞬時(shí)單位線與適當(dāng)?shù)漠a(chǎn)流模型結(jié)合，就能實(shí)現(xiàn)高板橋水庫(kù)流域的水庫(kù)洪水預(yù)報(bào)。

5 結(jié)論和討論

(1)本文以雁江區(qū)高板橋水庫(kù)流域?yàn)閷?shí)例，基于DEM提取出流域河網(wǎng)和Horton地貌參數(shù)，并成功建立了水庫(kù)流域的R-V地貌瞬時(shí)單位線，為進(jìn)一步探索小型水庫(kù)流域分布式水文預(yù)報(bào)模型奠定了基礎(chǔ)，對(duì)推動(dòng)小型水庫(kù)洪水預(yù)報(bào)方法研究具有積極意義。

(2)本研究采用中小流域暴雨洪水推理公式計(jì)算流域平均流速，計(jì)算公式中不僅包含下墊面特性參數(shù)，還充分考慮了降雨強(qiáng)度對(duì)流域平均流速的影響，且其參數(shù)值均可通過(guò)《四川省水文手冊(cè)》查算，具有一定的實(shí)用和推廣價(jià)值。

(3)本研究所使用的DEM為從互聯(lián)網(wǎng)免費(fèi)下載獲得，分辨率僅為30m，且經(jīng)過(guò)源數(shù)據(jù)插值處理，數(shù)據(jù)精度欠佳，以致提取出的流域參數(shù)與實(shí)際流域有一定出入。在實(shí)際的水庫(kù)洪水預(yù)報(bào)應(yīng)用中，為更準(zhǔn)確的獲取流域地形地貌信息，尚需采用精度更高的DEM產(chǎn)品。

(4)研究中發(fā)現(xiàn)，柵格河網(wǎng)生成時(shí)，設(shè)置不同的集水面積閾值，可得到不同密度的河網(wǎng)，主觀隨意性比較大?？蓢L試選取不同地區(qū)不同尺度的多個(gè)流域做專題研究，探尋河網(wǎng)生成閾值的固有特性，尋求不同區(qū)域不同尺度下的河網(wǎng)生成閾值常值。

(5)R-V GIUH理論未考慮河網(wǎng)的水動(dòng)力特性空間變異化，即認(rèn)為水質(zhì)點(diǎn)在流域中的傳播速度是固定不變的，然而實(shí)踐證明水流在河網(wǎng)中的傳播是一個(gè)非線性過(guò)程，且水流在從上游向下游傳播的過(guò)程中，其流速緩慢增加。同時(shí)，R-V GIUH理論幾乎完全忽略了坡面匯流過(guò)程，將坡面匯流時(shí)間忽略不計(jì)，這在坡面匯流滯時(shí)比重較大的小流域上，是無(wú)法真實(shí)模擬流域匯流過(guò)程的。故此，R-V地貌瞬時(shí)單位線用于小型水庫(kù)流域匯流計(jì)算，是存在一定理論缺陷的，尚需探索相關(guān)方面的改進(jìn)。