謝志剛，石 朋,2，紀(jì)小敏，趙蘭蘭，瞿思敏，崔彥萍

(1.河海大學(xué)水文水資源學(xué)院，南京 210098；2.河海大學(xué) 水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，南京 210098；3.江蘇省水文水資源勘測(cè)局，南京 210029；4.水利部信息中心，北京 100053)

已有研究表明[1,2]，流域匯流可被視為凈雨在地貌擴(kuò)散和水動(dòng)力擴(kuò)散綜合作用下的匯集過程。Redriguez-Iturbe和Gupta等[3,4]提出假定速度項(xiàng)表現(xiàn)為指數(shù)分布，則其線性反應(yīng)函數(shù)為時(shí)間指數(shù)分布，由此流域地貌瞬時(shí)單位線可以解釋為一滴雨抵達(dá)流域上任何處，其匯流時(shí)間的概率密度函數(shù)。芮孝芳等[5]研究發(fā)現(xiàn)，確定等待時(shí)間概率密度函數(shù)的實(shí)質(zhì)就是確定水質(zhì)點(diǎn)流速的概率密度函數(shù)，故把推求地貌瞬時(shí)單位線問題轉(zhuǎn)化為確定流速的概率密度函數(shù)問題。

計(jì)算流速的方法有很多，主要的流速計(jì)算公式包括曼寧公式、梅德門特公式、謝才公式等[6]。1993年，梅德門特教授[7]基于柵格假定了一個(gè)時(shí)不變的空間流速場(chǎng)；1995年，Muzik[8]在梅德門特教授的基礎(chǔ)上，結(jié)合曼寧公式與連續(xù)方程進(jìn)行流速計(jì)算，并利用GIS確定流速場(chǎng)；2003年，Assefa M.Melesse等[9]人總結(jié)出結(jié)合曼寧公式與運(yùn)動(dòng)波近似方程來估算坡面流速，聯(lián)合連續(xù)方程與曼寧公式求解河道流速，以此構(gòu)建流速場(chǎng)；2006年，石朋等[10,11]采用曼寧公式計(jì)算柵格單元流速并以此建立了一個(gè)在空間上變化而時(shí)間上不變的空間流速場(chǎng)，以沿渡河流域?yàn)檠芯繉?duì)象，取得了較好的計(jì)算效果?？追舱艿萚12]區(qū)分坡地單元與河道網(wǎng)格單元，以SCS公式計(jì)算坡地流速，以梅德門特公式計(jì)算河道流速，由此計(jì)算出整個(gè)流域的流速分布?？紤]到在一次降雨徑流過程中，雨強(qiáng)對(duì)流速的影響是不容忽視的[13]，本文擬通過將雨強(qiáng)納入考慮范圍并改進(jìn)原有的曼寧公式，計(jì)算流域各網(wǎng)格內(nèi)流速，并應(yīng)用于定安河流域，以提高洪水模擬過程的精度。

1 研究思路

1.1 流速計(jì)算方法

流域下墊面條件存在空間差異性，不同地形地貌條件下流域各處的流速大小與方向各不相同。流速大小的計(jì)算是本文的核心所在。本文主要選用了兩種流速計(jì)算方法，分別是：傳統(tǒng)曼寧公式及考慮雨強(qiáng)的曼寧公式。

1.1.1 曼寧公式

傳統(tǒng)的曼寧公式如下：

(1)

式中：n為地表粗糙系數(shù)，可以通過土地利用情況獲得；S為地表坡度，可以通過DEM獲得；R為水力半徑，可以近似以水深代替。

由已有研究可知[14]，水深的計(jì)算可以建立水深與匯水面積的函數(shù)關(guān)系：

H=φpAψp

(2)

式中：H為具有超越概率p的柵格平均水深；A為柵格上游流域匯水面積；φp為網(wǎng)格常數(shù)，與洪水頻率有關(guān)；ψp為幾何常數(shù)，與洪水頻率有關(guān)。

1.1.2 考慮雨強(qiáng)的曼寧公式

實(shí)際洪水過程中存在降雨過程，不同的降雨強(qiáng)度對(duì)匯流速度也存在影響，因此采用考慮雨強(qiáng)的曼寧公式進(jìn)行流速計(jì)算，可提高計(jì)算精度。計(jì)算公式如下：

(3)

式中：i為雨強(qiáng)大小的無量綱因子；b為經(jīng)驗(yàn)系數(shù)。

雨強(qiáng)作為影響流速大小的影響因子納入曼寧公式中，如公式(3)所示。在GIS中利用克里金插值法得出流域雨強(qiáng)柵格圖，即可知每個(gè)柵格內(nèi)的雨強(qiáng)大小，再由已知的糙率系數(shù)、坡度和水力半徑通過公式(3)可求出受雨強(qiáng)影響的柵格流速。由于整個(gè)洪水過程中存在降雨空間分布不均勻的情況[15]，不同地區(qū)雨強(qiáng)大小不同，并且在同一時(shí)段可能只有部分地區(qū)降雨，因此要對(duì)雨強(qiáng)柵格圖進(jìn)行0處理，即利用柵格計(jì)算器把所有0柵格替換成1，在無降雨地區(qū)i取值1，等價(jià)于直接用曼寧公式計(jì)算。

1.2 匯流時(shí)間計(jì)算

通過上述計(jì)算所得柵格流速，依據(jù)柵格大小可以計(jì)算出水流在柵格內(nèi)的滯留時(shí)間。計(jì)算公式如下：

Δτ=L/v

(4)

(5)

式中：Δτ為柵格內(nèi)滯留時(shí)間；L為柵格的邊長(zhǎng)。

公式(4)表示的是水流方向同柵格邊線平行的情況，公式(5)表示的是水流方向同柵格對(duì)角線平行的情況。在GIS中通過D8算法可以確定水流方向，即確定每個(gè)柵格內(nèi)水流方向，通過公式(4)和公式(5)計(jì)算各柵格內(nèi)滯留時(shí)間，在GIS中沿水流方向累積各柵格到達(dá)流域出口的時(shí)間即可確定各柵格匯流累積時(shí)間，計(jì)算公式如下：

(6)

式中：τ為某柵格到達(dá)流域出口的累積匯流時(shí)間；m為某柵格徑流路徑上的柵格數(shù)。

1.3 單位線的確定

通過以上兩種流速計(jì)算方法在GIS中可利用柵格計(jì)算器計(jì)算出每個(gè)柵格的流速大小，以此可得出兩種不同的空間流速場(chǎng)，利用空間流速場(chǎng)以公式(4)、(5)可以計(jì)算得到每個(gè)柵格的滯留時(shí)間，通過公式(6)沿水流方向得出各柵格匯流累積時(shí)間，在GIS中利用重分類模塊統(tǒng)計(jì)時(shí)段內(nèi)通過流域出口的柵格數(shù)，即得到時(shí)段內(nèi)出流面積。以匯流時(shí)間為橫坐標(biāo)，以時(shí)段內(nèi)出流面積與總面積比值為縱坐標(biāo)，可得到流域匯流時(shí)間的概率密度分布。根據(jù)徑流過程形成的“粒子學(xué)說”[2-4]，即水質(zhì)點(diǎn)在弱相互作用情況下的匯流時(shí)間概率密度分布函數(shù)等價(jià)于地貌瞬時(shí)單位線，可知以上所得匯流時(shí)間概率密度分布即為地貌瞬時(shí)單位線，將此應(yīng)用于實(shí)例當(dāng)中，對(duì)比徑流模擬計(jì)算結(jié)果，分析兩種流速計(jì)算方法存在的問題與優(yōu)點(diǎn)。

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

2.1 流域資料處理

本文所選流域是萬泉河水系內(nèi)的定安河流域，流域面積1 333 km2，內(nèi)設(shè)長(zhǎng)征、大墩、合羅、加報(bào)、羅擔(dān)、木色、瓊中、石古、思河、烏坡和烏石11個(gè)雨量站，具有完善的降雨資料。流域地處熱帶季風(fēng)氣候區(qū)，常年雨量充沛，多年平均降雨量約1 639 mm。流域是萬泉河水系的一級(jí)支流，受水利工程影響較小，具有多年連續(xù)的降雨徑流資料，適合本文模擬研究計(jì)算。

本文基于流域DEM(90×90)進(jìn)行空間流速場(chǎng)的構(gòu)建，推導(dǎo)出相應(yīng)的地貌瞬時(shí)單位線，并應(yīng)用于降雨徑流模擬分析。流域坡度分布如圖1所示。

圖1 流域坡度分布Fig.1 Slope distribution of River Basin

根據(jù)流域土地利用情況確定了糙率系數(shù)n的柵格數(shù)據(jù)圖；利用GIS中水文分析模塊可計(jì)算各柵格上游匯水面積，以此根據(jù)公式(2)可求出流域水深柵格數(shù)據(jù)圖；通過克里金插值法得出流域時(shí)段雨強(qiáng)柵格數(shù)據(jù)圖。利用GIS中柵格計(jì)算器，以上述得出的柵格數(shù)據(jù)通過公式(1)～(3)計(jì)算每個(gè)柵格的流速，分別構(gòu)建出以曼寧公式計(jì)算出的空間流速場(chǎng)和以考慮雨強(qiáng)的曼寧公式計(jì)算出的空間流速場(chǎng)，雨強(qiáng)曼寧公式以20010825次洪為例。圖2(a)為曼寧公式構(gòu)建出空間上變化而時(shí)不變流速場(chǎng)，圖2(b)為2001年8月25日次洪降雨第2個(gè)時(shí)段雨強(qiáng)的空間流速場(chǎng)柵格圖。

圖2 計(jì)算所得空間流速場(chǎng)Fig.2 Calculated spatial velocity field

根據(jù)匯流時(shí)間計(jì)算公式(4)～(6)可以計(jì)算出相應(yīng)的匯流累積時(shí)間柵格分布圖，圖3(a)為曼寧公式計(jì)算出的累積時(shí)間柵格圖，圖3(b)為2001年8月25日次洪降雨第2個(gè)時(shí)段雨強(qiáng)的計(jì)算出的累積時(shí)間柵格圖。

圖3為沿水流方向到達(dá)流域出口的匯流累積時(shí)間分布圖，距離匯流出口越遠(yuǎn)所需匯流時(shí)間越長(zhǎng)，與實(shí)際情況相符合，據(jù)此統(tǒng)計(jì)匯流時(shí)間的概率密度分布函數(shù)[16]，得出地貌瞬時(shí)單位線。雨強(qiáng)曼寧公式由于各時(shí)段雨強(qiáng)不同所得出的單位線也不同，故采用時(shí)段出流方式進(jìn)行匯流計(jì)算。

圖3 匯流累積時(shí)間分布Fig.3 Accumulation time distribution of confluence

2.2 降雨徑流模擬計(jì)算

利用上述所得地貌瞬時(shí)單位線轉(zhuǎn)換成1 h時(shí)段單位線，應(yīng)用于定安河流域進(jìn)行降雨徑流模擬分析，選取該流域2000-2013年中的7場(chǎng)洪水進(jìn)行模擬計(jì)算，流域每年洪澇災(zāi)害多發(fā)生于5-11月之間，強(qiáng)降雨主要集中在7-10月，考慮雨強(qiáng)對(duì)流速的影響，故選取此期間7場(chǎng)洪水進(jìn)行模擬計(jì)算，所得計(jì)算結(jié)果如表1所示。經(jīng)過對(duì)比分析，由曼寧公式計(jì)算值相比較于考慮雨強(qiáng)的曼寧公式計(jì)算值，7場(chǎng)洪水平均徑流相對(duì)誤差由11.52%下降到5.93%，平均洪峰相對(duì)誤差由6.32%下降到3.85%，平均確定性系數(shù)由0.761上升到0.840。由此可知考慮雨強(qiáng)大小影響的曼寧公式計(jì)算所得結(jié)果更符合實(shí)際情況，其徑流深、洪峰流量、確定性系數(shù)模擬效果更優(yōu)[17]。圖4為20010825和20071011兩場(chǎng)洪水由曼寧公式計(jì)算的洪水過程線和考慮雨強(qiáng)的曼寧公式計(jì)算的洪水過程性與實(shí)測(cè)洪水過程線的對(duì)比圖。

表1 定安河徑流計(jì)算結(jié)果Tab.1 Runoff calculation results of the Ding An River

注：此平均值為絕對(duì)值平均值。

圖4 洪水模擬計(jì)算對(duì)比Fig.4 Comparison of flood simulation

3 結(jié) 語

本文核心是在流域流速計(jì)算方法上，而傳統(tǒng)地貌瞬時(shí)單位線的確定是建立在流速空間分布均勻的情況下，實(shí)際流域中流速大小受地形地貌影響以及水力條件影響較大，故本文提出了考慮雨強(qiáng)影響的曼寧公式計(jì)算流域流速，建立了整個(gè)流域的空間流速場(chǎng)。曼寧公式充分考慮了土地利用情況、坡度以及水深對(duì)流速的影響，由此可以確定一個(gè)空間上變化而時(shí)不變的空間流速場(chǎng)，但是實(shí)際洪水過程中存在降雨空間分布不均勻性，不同區(qū)域雨強(qiáng)不同亦對(duì)流速產(chǎn)生重大影響，故此加入雨強(qiáng)這一變量因素，以此改進(jìn)曼寧公式求得流域空間流速場(chǎng)，應(yīng)用于定安河流域，取得了比較好的模擬效果，故此可推廣應(yīng)用于缺資料地區(qū)的匯流計(jì)算過程中。

本文所取雨強(qiáng)為1 h時(shí)段雨強(qiáng)，未考慮不同時(shí)段大小雨強(qiáng)對(duì)流速的影響，后期重點(diǎn)研究不同時(shí)段大小雨強(qiáng)對(duì)流速計(jì)算精度的影響，提高徑流模擬效果。本文所取柵格大小是90 m×90 m，未考慮不同分辨率大小情況下坡度值精度不同，由此產(chǎn)生流速計(jì)算誤差，故在后續(xù)研究過程中會(huì)進(jìn)一步考慮分辨率對(duì)流速計(jì)算產(chǎn)生的影響。