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

包紅軍，王莉莉，李致家

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

0 引 言

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

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

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

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

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

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

K=ΔL/VW

(1)

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

(2)

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

(3)

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

Q0=Qm+0.5(QP-Qm)

(4)

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

對(duì)于給定流量過(guò)程可以用公式(4)來(lái)推求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等于過(guò)水?dāng)嗝娴乃鎸扺。

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

A=2yW/3

(6)

式中，y為水深。對(duì)于拋物線型斷面水力半徑可以根據(jù)Koegelenberg在1997年研究成果[21](見(jiàn)表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ù)河段長(zhǎng)度由式(1)確定K的值。

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

(12)

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

河道洪水預(yù)報(bào)包括流量演算法和水位計(jì)算法。行蓄洪是通過(guò)流量的變化(分洪、蓄洪)來(lái)影響水位的。因而只要在流量演算時(shí)考慮行蓄洪的影響即可。以河段[i，i+1]為例，見(jiàn)圖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]的來(lá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為變化后的水位(見(jiàn)圖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)相對(duì)較小[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è)主要支流匯入，另外還有洪河分洪道、谷河和潤(rùn)河匯入和蒙洼蓄洪區(qū)、姜唐聯(lián)湖蓄洪區(qū)、南潤(rùn)段、邱家湖、潤(rùn)趙段三個(gè)行洪區(qū)。王家壩至魯臺(tái)子河段共有三個(gè)水文站與四個(gè)水位站，在整個(gè)流域防汛中，一直處于洪水預(yù)報(bào)關(guān)鍵之處[24-26]。

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

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

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

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

3 結(jié)論與討論

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

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

(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)誤差，建立對(duì)系統(tǒng)模型與預(yù)報(bào)的現(xiàn)時(shí)校正的回饋機(jī)制，是提高精度的重要手段[7]。

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