基于Muskingum-Cunge法的河道水位流量預(yù)報(bào)研究

包紅軍，王莉莉，李致家

(1.國家氣象中心，北京100081；2.河海大學(xué)水文水資源學(xué)院， 江蘇南京210098)

0 引 言

河道水流演算水力要素屬于三維非恒定的水力學(xué)問題[1]。考慮到三維非恒定的水力學(xué)模型基本方程理論假設(shè)與數(shù)學(xué)求解問題，在生產(chǎn)實(shí)際作業(yè)中往往概化為一維非恒定流模型[2- 6]。圣維南方程組是描述一維水流非恒定流運(yùn)動(dòng)的基礎(chǔ)方程，水位模擬與預(yù)報(bào)精度較高；但其對資料，特別是斷面資料要求較高，計(jì)算繁瑣，且目前尚無法求出解析通解，其近似概化模型有運(yùn)動(dòng)波、擴(kuò)散波等[7]。其他用于河道洪水演算方法還有滯后演算法等經(jīng)驗(yàn)方法、線性水庫法、非線性水庫法和Muskingum法，其中以Muskingum最為常用[8]。

Muskingum法是實(shí)際生產(chǎn)中常用的集總式河道演算方法，由于在河道洪水演算中的簡便性和廣泛的適用性，其參數(shù)確定方法也多種多樣。最初的Muskingum法參數(shù)槽蓄系數(shù)K和權(quán)重因子x的確定方法是試錯(cuò)法[1-2]。為了克服試錯(cuò)法計(jì)算量大和主觀性過強(qiáng)的缺點(diǎn)，有學(xué)者采用最小二乘法推導(dǎo)出參數(shù)K和x的計(jì)算公式[9]。另外，Singh對Muskingum法參數(shù)確定的方法歸納為：最小二乘法或者圖解法、矩陣及累積量法、直接試解優(yōu)選法[10-11]。最優(yōu)化方法實(shí)質(zhì)上屬于“黑箱子”模型，難以保證K和x的物理意義[12-13]。1969年，Cunge證明了Muskingum法是擴(kuò)散波的二階近似解，將Muskingum法與水力學(xué)中的擴(kuò)散波聯(lián)系起來，提出Muskingum-Cunge法；另一方面，其形式與Muskingum法相同，但有明確的水力學(xué)意義[14-15]。

本研究以淮河中游河道為例，對Muskingum-Cunge法參數(shù)進(jìn)行基于河道流量演算中的經(jīng)驗(yàn)關(guān)系推求，并且結(jié)合水文水位法、擴(kuò)散波非線性水位法，改善Muskingum-Cunge法不能模擬水位的問題，形成基于Muskingum-Cunge法的河道水位流量預(yù)報(bào)模型，在淮河干流河道洪水預(yù)報(bào)中進(jìn)行驗(yàn)證，以探討減少方法對資料的依賴程度。

1 基于Muskingum-Cunge法的河道水位流量預(yù)報(bào)模型

1.1 Muskingum-Cunge參數(shù)推求

Muskingum-Cunge法中參數(shù)K表示洪水波在河段長為ΔL中的傳播運(yùn)動(dòng)歷時(shí)，假定洪水波運(yùn)動(dòng)速度為VW，則K可由下式計(jì)算[16-17]。即

K=ΔL/VW

(1)

1989年，Wilson和Ruffini研究得到：對于某個(gè)特定矩形、三角形和拋物線河道斷面，其波速VW分別是斷面平均波速Vav的5/3、4/3和11/9[17]。而根據(jù)曼寧公式斷面平均流速[18]

(2)

式中，n為曼寧糙率系數(shù)，無因次；R為水力半徑，m；S為水面比降，m/m。由此可見，當(dāng)推求出R值，即可得到斷面平均流速Vav。式(2)可寫成[19]

(3)

式中，A為斷面的過水面積，m2；Q0為參考流量，m3/s(Reference Discharge)。Q0值求解根據(jù)文獻(xiàn)[17，19]，采用式(4)來求解某場洪水中的Q0。即

Q0=Qm+0.5(QP-Qm)

(4)

式中，Qm為過程最小流量，m3/s；QP為過程洪峰流量，m3/s。

對于給定流量過程可以用公式(4)來推求Q0。糙率n和坡度S分別取自河道地形資料。Punmia和Pande[20]在1981年研究得出在自然流域中濕周與流量的經(jīng)驗(yàn)關(guān)系。即

(5)

式中，P為濕周；c為系數(shù)，值域在4.71～4.81。寬淺河道時(shí)，可認(rèn)為P等于過水?dāng)嗝娴乃鎸扺。

拋物線型斷面面積可近似按下述公式計(jì)算[18]

A=2yW/3

(6)

式中，y為水深。對于拋物線型斷面水力半徑可以根據(jù)Koegelenberg在1997年研究成果[21](見表1)計(jì)算

R?d=2y/3

(7)

式中，d為平均水深，m。

表1 水力半徑與水深經(jīng)驗(yàn)關(guān)系

將式(5)、(6)代入式(3)，可得

(8)

因?yàn)闈裰躊在寬淺河道中可近似為水面寬W，所以拋物線型斷面的水深

(9)

同理可得到當(dāng)斷面形狀為矩形時(shí)的水深

(10)

以及當(dāng)斷面形狀為三角形時(shí)水深

(11)

根據(jù)表1可得到各種河道斷面形狀的水力半徑R，代入式(11)求得Vav；同時(shí)，也得出Vav是坡度S單調(diào)增函數(shù)，隨著S的增加而增加的結(jié)果；進(jìn)而，根據(jù)河段長度由式(1)確定K的值。

Muskingum-Cunge法參數(shù)x根據(jù)1993年Fread[16]提出的公式計(jì)算。即

(12)

1.2 考慮行蓄洪區(qū)的Muskingum-Cunge法

河道洪水預(yù)報(bào)包括流量演算法和水位計(jì)算法。行蓄洪是通過流量的變化(分洪、蓄洪)來影響水位的。因而只要在流量演算時(shí)考慮行蓄洪的影響即可。以河段[i，i+1]為例，見圖1。

圖1 河段[i，i+1]概化示意

Muskingum-Cunge法遞推公式為

(13)

根據(jù)圖1，考慮行蓄洪區(qū)時(shí)，i斷面入流包括4個(gè)：①支流匯入Gi；②行洪分流出入流(分左右兩側(cè))QFli，QFfi；③蓄洪入流(分左右兩側(cè))QSli、QSfi；(4)[i-1，i]的來水量；則

(14)

(15)

考慮行蓄洪區(qū)的Muskingum-Cunge法遞推公式

(16)

行洪區(qū)與蓄洪區(qū)處理方式采用文獻(xiàn)[22]的做法：行洪區(qū)進(jìn)行一維河道洪水演算，蓄洪區(qū)一般有閘控制，不進(jìn)行水流演算。

1.3 基于柱蓄和楔蓄的水文水位法

設(shè)定1-1為初始水位，2-2為變化后的水位(見圖2)，則河道槽蓄量增量

dW=BLdZ下+BLX1(dZ上-dZ下)

(17)

(18)

圖2 柱蓄和楔蓄示意

(19)

式中，i0為河道河底比降；iΔ為水流附加比降。

(20)

與Muskingum-Cunge法聯(lián)解，實(shí)現(xiàn)斷面水位推求如下

(21)

(22)

(23)

(24)

(25)

1.4 擴(kuò)散波非線性水位法

利用基于柱蓄和楔蓄的水文水位法與Muskingum-Cunge法聯(lián)解求得逐段流量和邊界水位后，可采用圣維南方程組中的動(dòng)量方程求解水位[22]。在淮河干流，慣性項(xiàng)相對較小[7]，可忽略去慣性項(xiàng)，得到

(26)

(27)

以θ為差分因子，應(yīng)用四點(diǎn)隱式差分格式，得

(28)

再利用迭代計(jì)算推求河道斷面水位。

2 基于Muskingum-Cunge法的河道水位流量預(yù)報(bào)模型應(yīng)用分析

本次研究以王家壩水文站至魯臺(tái)子水文站的淮河干流為試驗(yàn)河段。河段右岸有史河、淠河兩個(gè)支流匯入，并且還有兩個(gè)蓄洪區(qū)(城西湖蓄洪區(qū)、城東湖蓄洪區(qū))；左岸有潁河一個(gè)主要支流匯入，另外還有洪河分洪道、谷河和潤河匯入和蒙洼蓄洪區(qū)、姜唐聯(lián)湖蓄洪區(qū)、南潤段、邱家湖、潤趙段三個(gè)行洪區(qū)。王家壩至魯臺(tái)子河段共有三個(gè)水文站與四個(gè)水位站，在整個(gè)流域防汛中，一直處于洪水預(yù)報(bào)關(guān)鍵之處[24-26]。

王家壩至魯臺(tái)子河段全長155.16 km。根據(jù)王家壩～魯臺(tái)子河道水利工程的位置、支流的匯入位置、水文站水位站的位置和行蓄洪區(qū)位置情況，河段共分10段。

選取淮河1996年至2008年汛期洪水資料進(jìn)行檢驗(yàn)。14場洪水均取得較好的模擬效果。為了能使Musking-Cunge能夠模擬好大洪水時(shí)干流洪水頂托作用，支流延后至退水期再進(jìn)入干流的洪水過程，故假設(shè)在支流洪水進(jìn)入干流之前，存在一個(gè)線性水庫，即支流洪水需要經(jīng)過一虛擬線性水庫調(diào)蓄后方能進(jìn)入干流。對于綜合法而言，當(dāng)河道比降大時(shí)，屬于運(yùn)動(dòng)波，用水文水位法提供上邊界水位條件，而當(dāng)河底比降小時(shí)，屬于擴(kuò)散波，用水文水位法提供下邊界條件。前者用擴(kuò)散波非線性水位法從上游向下游推求其他斷面的水位，后者用擴(kuò)散波非線性水位法從下游向上游推求其他斷面的水位。淮河干流河道比降較低，因而采用擴(kuò)散波非線性水位法自下游向上游推求水位。通過在1996年至2008年洪水的檢驗(yàn)，該處理方式具有一定的效果，使得洪水的模擬精度得到了提高，如表2所示。

表2 王家壩至魯臺(tái)子河段洪水模擬結(jié)果

注：y表示行蓄洪區(qū)使用年份；①為實(shí)測最高水位減去計(jì)算最高水位；②為計(jì)算與實(shí)測水位過程擬合的確定性系數(shù)；③為計(jì)算洪峰流量的相對誤差；④為計(jì)算與實(shí)測流量過程擬合的確定性系數(shù)。

3 結(jié)論與討論

本次研究對Muskingum-Cunge法參數(shù)進(jìn)行基于河道流量演算中的經(jīng)驗(yàn)關(guān)系推求，并且結(jié)合基于柱蓄與楔蓄的水文水位法、擴(kuò)散波非線性水位法，改善Muskingum-Cunge法不能模擬水位的問題，建立基于Muskingum-Cunge法的河道水位流量預(yù)報(bào)模型。通過在淮河干流河道洪水模擬中驗(yàn)證，效果較好。

(1)Muskingum是根據(jù)實(shí)測洪水資料優(yōu)選參數(shù)，行蓄洪等水利工程以及河道洪水自身特性等對河道洪水的影響可以在優(yōu)選出來的參數(shù)值中反映出來，而Muskingum-Cunge法的參數(shù)完全根據(jù)水力要素推求，沒有考慮到分洪、漫灘等等，推求出來的參數(shù)不能反映水利工程和河道自身水流特性等等因素的影響。這是Muskingum-Cunge法在下一步研究與實(shí)際生產(chǎn)應(yīng)用中需要解決的問題之一。

(2)為了實(shí)現(xiàn)Muskingum-Cunge法模擬水位，Muskingum-Cunge法與基于柱蓄與楔蓄的水文水位法、擴(kuò)散波非線性水位形成基于Muskingum-Cunge法的河道水位流量預(yù)報(bào)模型，方法的優(yōu)點(diǎn)在于不依賴洪水歷史資料，參數(shù)完全基于河道物理特征求得。為了提高水位模擬精度，應(yīng)用實(shí)時(shí)校正模型，利用當(dāng)前的預(yù)報(bào)誤差，建立對系統(tǒng)模型與預(yù)報(bào)的現(xiàn)時(shí)校正的回饋機(jī)制，是提高精度的重要手段[7]。

對于復(fù)雜的分叉水系及分紅、潰口、潰壩等的洪水預(yù)報(bào)，或者在計(jì)算河道內(nèi)興建了水利工程或者河道特征發(fā)生了明顯變化后，需要計(jì)算和預(yù)報(bào)河道內(nèi)任何斷面任何時(shí)刻的水位、流速、流量等要素時(shí)，這時(shí)水動(dòng)力學(xué)方法仍是首選[6，25]。