基于Muskingum-Cunge法的河道水位流量預報研究

包紅軍，王莉莉，李致家

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

0 引 言

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

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

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

1 基于Muskingum-Cunge法的河道水位流量預報模型

1.1 Muskingum-Cunge參數推求

Muskingum-Cunge法中參數K表示洪水波在河段長為ΔL中的傳播運動歷時，假定洪水波運動速度為VW，則K可由下式計算[16-17]。即

K=ΔL/VW

(1)

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

(2)

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

(3)

式中，A為斷面的過水面積，m2；Q0為參考流量，m3/s(Reference Discharge)。Q0值求解根據文獻[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年研究得出在自然流域中濕周與流量的經驗關系。即

(5)

式中，P為濕周；c為系數，值域在4.71～4.81。寬淺河道時，可認為P等于過水斷面的水面寬W。

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

A=2yW/3

(6)

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

R?d=2y/3

(7)

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

表1 水力半徑與水深經驗關系

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

(8)

因為濕周P在寬淺河道中可近似為水面寬W，所以拋物線型斷面的水深

(9)

同理可得到當斷面形狀為矩形時的水深

(10)

以及當斷面形狀為三角形時水深

(11)

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

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

(12)

1.2 考慮行蓄洪區的Muskingum-Cunge法

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

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

Muskingum-Cunge法遞推公式為

(13)

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

(14)

(15)

考慮行蓄洪區的Muskingum-Cunge法遞推公式

(16)

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

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

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

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

(17)

(18)

圖2 柱蓄和楔蓄示意

(19)

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

(20)

與Muskingum-Cunge法聯解，實現斷面水位推求如下

(21)

(22)

(23)

(24)

(25)

1.4 擴散波非線性水位法

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

(26)

(27)

以θ為差分因子，應用四點隱式差分格式，得

(28)

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

2 基于Muskingum-Cunge法的河道水位流量預報模型應用分析

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

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

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

表2 王家壩至魯臺子河段洪水模擬結果

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

3 結論與討論

本次研究對Muskingum-Cunge法參數進行基于河道流量演算中的經驗關系推求，并且結合基于柱蓄與楔蓄的水文水位法、擴散波非線性水位法，改善Muskingum-Cunge法不能模擬水位的問題，建立基于Muskingum-Cunge法的河道水位流量預報模型。通過在淮河干流河道洪水模擬中驗證，效果較好。

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

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

對于復雜的分叉水系及分紅、潰口、潰壩等的洪水預報，或者在計算河道內興建了水利工程或者河道特征發生了明顯變化后，需要計算和預報河道內任何斷面任何時刻的水位、流速、流量等要素時，這時水動力學方法仍是首選[6，25]。