劉力真,李福仲，陸小蕾

（1.山東省水利科學(xué)研究院,山東 濟南 250013；2.山東大學(xué)土建與水利學(xué)院，山東 濟南 250061）

目前，流域匯流單位線推求多是通過兩種方法：一是依據(jù)現(xiàn)有水文資料推求；二是依據(jù)典型單位線結(jié)合流域?qū)嶋H情況進(jìn)行縮放。如此，對反應(yīng)實際匯流情況均有一定偏差；前者過多依賴于水文資料，對資料的準(zhǔn)確性和長序列要求過高，后者回避了流域單位線的多樣性，與實際不符?；诘乩硇畔①Y料，推求符合流域地理特性的單位線，解決由于暴雨中心不同引起單位線變化的問題。

1 原 理

1.1 流域出口斷面流量的形成

一場降雨在一個網(wǎng)格中形成的產(chǎn)流量可視為一個“水滴”，水滴的大小可通過調(diào)整網(wǎng)格的尺寸來設(shè)定，這樣定義的“水滴”可稱為“網(wǎng)格水滴”。并且“水滴”降落在流域各處是隨機的，扣除損失后將會選擇一條流路、花費一定時間最終流達(dá)流域出口斷面，其所花費的時間即為該水滴的流域匯流時間。

假定流域匯流系統(tǒng)是線性系統(tǒng)，忽略相鄰網(wǎng)格水流之間的水力聯(lián)系，即各個網(wǎng)格水滴分別匯流，互不干擾，只在流域出口處做同時到達(dá)的水滴量累加即得到流域出口的流量過程。

假設(shè)一場降雨由n個水滴組成，其編號分別為1，2，…，i，…n。對其中任一水滴 i，若其流達(dá)流域出口斷面的流路長為ri，通過該流路的運動速度一般與流路有關(guān)，記為vi(x)，則該水滴的流域匯流時間為：

式中：τi為第i個水滴的流域匯流時間；x為流路。

對每個水滴使用式(1)，就可得到水滴流域匯流時間的空間分布。若降雨是在t=0時刻降落到流域上的，則不難理解，只有那些流域匯流時間正好等于T的水滴才能在t=T時刻到達(dá)流域出口斷面，組成T時刻的流域出口斷面流量，所以有

式中：q(t)為t=0時刻瞬時降雨形成的t=T時刻流域出口斷面流量；

I(τi=t)為指示函數(shù)，當(dāng) τi=t時，其值為 1，表明屬于τi=t的水滴是組成T時刻出口斷面流量的一分子，否則為0，表明不屬于τi=t的水滴不能成為T時刻出口斷面流量的一分子；

ΔAi為第i個水滴占據(jù)的面積；

hi為t=0時刻的降雨強度。

對于一次降雨過程，則可通過在不同降雨時刻重復(fù)使用式(2)，然后利用疊加原理來推求出口斷面流量過程線。

應(yīng)用上述方法推求降雨形成的流域出口斷面流量過程的關(guān)鍵是如何確定水滴尺度及其流路和速度。

1.2 D8法確定水滴尺度及流路

一場降雨在一個網(wǎng)格中形成的產(chǎn)流量可視為一個“水滴”，水滴的大小可通過調(diào)整網(wǎng)格的尺寸來設(shè)定，這樣定義的“水滴”可稱為“網(wǎng)格水滴”。

每個網(wǎng)格水滴的流路應(yīng)用數(shù)字高程模型（DEM）進(jìn)行流向分析。如圖1所示，在網(wǎng)格化的流域中，與其中任一網(wǎng)格相鄰的網(wǎng)格共有8個，可認(rèn)為任一網(wǎng)格中的產(chǎn)流量向相鄰網(wǎng)格流動有8種可能。假設(shè)其中最大坡度方向就是實際發(fā)生的水流方向，只要計算出一個網(wǎng)格水滴的8個方向，則其中最大者就是水流方向。這種確定流向的方法稱為D8法。

圖1 D8法示意圖

任一網(wǎng)格與其相鄰8個網(wǎng)格之間的坡度可按下式計算：

式中：Sij為第i個網(wǎng)格與相鄰網(wǎng)格中任一個之間的坡度；Zj為第i個網(wǎng)格中心之高程；

Zj為與第i個網(wǎng)格相鄰的第j個網(wǎng)格中心之高程；

d為網(wǎng)格的邊長。

流域上每個網(wǎng)格水滴的流向確定后，就可求得每個網(wǎng)格水滴流路的走向及其幾何長度。

上述D8法屬于單流向法。此外還有多流向法可供選用。經(jīng)驗表明，對于流域匯流計算采用D8法已能獲得令人滿意的精度。

1.3 流速的確定

1）坡面流速公式

式中：SS為坡面坡度；a為反映坡面糙率的經(jīng)驗系數(shù)。

2）河道流速公式

式中：A為網(wǎng)格以上匯水面積；b為經(jīng)驗指數(shù)；其余符號的意義同前述。

1.4 網(wǎng)格匯流時間計算

每一個網(wǎng)格單元的流速經(jīng)式(4)或式(5)計算后,可由式(6)計算徑流在每一個網(wǎng)格內(nèi)的滯時,即

式中：L為網(wǎng)格的邊長。

若沿著徑流路徑,“水滴”經(jīng)過n個網(wǎng)格達(dá)到出口，可由式(7)計算出每一個網(wǎng)格單元到達(dá)流域出口的匯流時間,即

2 引入暴雨中心矩陣

考慮到暴雨中心常處在變化中，導(dǎo)致雨量分布不均勻，匯流時間也有所不同。為提高匯流計算的精度，本次算法引入了暴雨中心矩陣：

1）矩陣元素平均數(shù)為1；

2）暴雨中心位置的矩陣元素數(shù)值最大，離暴雨中心最遠(yuǎn)的流域位置矩陣元素數(shù)值最??；

3）矩陣元素的數(shù)值與離暴雨中心的距離成反比。

3 地貌單位線在沿渡河流域的應(yīng)用

3.1 DEM預(yù)處理

將引入暴雨矩陣的瞬時單位線計算方法應(yīng)用到沿渡河流域。以長江流域內(nèi)的沿渡河流域為研究對象。利用GIS工具對DEM填洼后,提取流域流向、邊界、坡度、網(wǎng)格集水面積和水系等，并以矩陣形式存儲。計算流域中每一個網(wǎng)格的坡度和匯水面積，同時對坡度為0的網(wǎng)格進(jìn)行處理。采取Matlab軟件編程運用D8法可以得到任一網(wǎng)格沿匯流路徑到達(dá)流域出口的距離。

3.2 參數(shù)率定

選擇流速公式（4）、（5），給定其參數(shù)初值，計算流域匯流時間矩陣，并統(tǒng)計其概率分布，得地貌瞬時單位線。結(jié)合地面凈雨計算結(jié)果，進(jìn)行流域匯流模擬，以計算流量過程線和實測流量過程線擬合最佳為目標(biāo)函數(shù)，率定參數(shù)a和b。

采用Matlab編程計算，共模擬4場洪水過程，分別得出每場洪水的a、b值，采用算術(shù)平均法得出最終結(jié)果（見表 1）。

表1 沿渡河流域洪水調(diào)參過程參數(shù)表

調(diào)參結(jié)果：a=0.954，b=0.9

推算暴雨中心矩陣的計算思路：

其中 c為未知系數(shù)，d為zx（i，j）距離暴雨中心的距離，i和j分別為流域的橫向和縱向元素個數(shù)。

4 結(jié)語

基于“網(wǎng)格水滴”的流域匯流計算方法，是一個理論依據(jù)較強的一般性流域匯流計算方法。既考慮了流域形狀、水系分布等對流域匯流的影響，又可通過坡度分布考慮空間流速場對流域匯流的作用；既保留了面積-時間曲線的優(yōu)點，又克服了單位線以降雨空間分布均勻為前提的缺點，既可用于分布式流域匯流，也可用于集總式流域匯流。

[1] 任立良,劉新仁.數(shù)字高程模型信息提取與數(shù)字水文模型研究進(jìn)展[J].水科學(xué)進(jìn)展[J],2000,11(4):463-469.

[2] 李紀(jì)人.地理信息系統(tǒng)在水利中的應(yīng)用[J].中國水利,2001(8):167-168.

[3] 左其亭,王中根.現(xiàn)代水文學(xué)[M].鄭州:黃河水利出版社,2002.

[4] 王中根,劉昌明,吳險峰.基于DEM的分布式水文模型研究綜述[J].自然資源學(xué)報，2003,18(2):168-173.