河道匯流演算方法對比分析與研究

龔 定，李致家，臧帥宏，孫明坤

(河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098)

河道匯流演算實(shí)質(zhì)上是一個水力學(xué)問題，由于存在需要較多的河道實(shí)測斷面資料與較長的計算耗時等問題，水文學(xué)家發(fā)展了水文學(xué)方法和簡化的水力學(xué)方法[1-3]。1938年McCarthy提出一種原始的水文洪水演算方法，即Muskingum法[4-5]，將所有的地形特征和河段的水力學(xué)特性歸為模型的系列線性參數(shù)[6]。20世紀(jì)60年代，趙人俊[1]推導(dǎo)出分段馬斯京根法連續(xù)演算的通用公式，解決了長河道預(yù)見期短的問題。實(shí)際河道水流條件復(fù)雜多變，河道的水力特性可能隨沿程發(fā)生變化，又形成了各種非線性的演算方法。1969年Cunge在Muskingum法的基礎(chǔ)上提出了Muskingum-Cunge(MC)法[7]；1978年P(guān)ouce等在Cunge法的基礎(chǔ)上，考慮波速的變動問題，將Muskingum演算法擴(kuò)展到時間可變參數(shù)[8]，提出了變動參數(shù)MC演算法(以下簡稱MC演算法)。使用線性參數(shù)的Muskingum法很好地保證了演算過程中河段的蓄量平衡，但MC演算法演算過程導(dǎo)致的質(zhì)量損失隨著河底坡度的減小而增加[9-10]。針對MC演算法中的質(zhì)量不守恒問題，2007年Todini等[11-14]提出Muskingum-Cunge-Todini(MCT)可變參數(shù)法，引入了修正因子和校正系數(shù)，進(jìn)一步完善了河道匯流演算方法。

河道匯流演算在流域洪水預(yù)報中起著關(guān)鍵性的作用。我國幅員遼闊，流域特點(diǎn)各不相同，匯流條件復(fù)雜[15-16]。在匯流過程中，河道內(nèi)水流運(yùn)動受地形影響較大[17]，傳統(tǒng)的河道演算方法難以準(zhǔn)確地描述河道的匯流特性，降低了洪水預(yù)報的精度[18-20]。MC演算法和MCT可變參數(shù)法在國外已經(jīng)取得較為廣泛的應(yīng)用[21-22]，在國內(nèi)，非線性和變參數(shù)河道匯流演算法在生產(chǎn)上還沒有較為成熟的理論和應(yīng)用[1]。

本文在分段馬斯京根法的基礎(chǔ)上選用考慮河道斷面形狀和水力特性的MC演算法和MCT可變參數(shù)法對沅水流域沅陵站—王家河站河段進(jìn)行河道匯流演算，對比分析考慮真實(shí)河道特性的匯流演算方法模擬效果。

1 方 法 介 紹

1.1 MC演算法

Ponce和Yevjevich在沒有改變問題實(shí)質(zhì)的基礎(chǔ)上[8]，提出了出流系數(shù)的表達(dá)式：

(1)

式中：C1、C2、C3為演算系數(shù)；C、D分別為庫朗數(shù)和單元格雷諾數(shù)，即物理和數(shù)字?jǐn)U散率；Δx為河流或河道被分割計算區(qū)段的長度，m；Δt為積分時間步長，s；B為地表寬度，m；S0為底坡坡度；Q為參比流量，m3/s；c為波速，m/s。

1.2 MCT可變參數(shù)法

MCT可變參數(shù)法并沒有改變MC演算法的簡單性，但同時又解決了物質(zhì)守恒的問題[23]。由式(1)可得

(2)

MC演算法中，在一個時間步長內(nèi)，k和ε不隨時間變化[24]。即在計算相鄰時刻的河段蓄量時，使用的k和ε是不一致的，在t時刻，有：

(3)

式中下標(biāo)代表相鄰兩個河段?？梢缘玫剑?/p>

(4)

式中St為t時刻的河段槽蓄量，上標(biāo)代表同一個河段蓄量的不同計算方式，這就造成了在連續(xù)時段下，蓄量出現(xiàn)不守恒的情況。在MCT可變參數(shù)法中，k和ε在一個時間步長內(nèi)變化，在每個時間步長分別產(chǎn)生兩個庫朗數(shù)和單元格雷諾數(shù)，解決了式(4)蓄量不平衡的問題。

MCT可變參數(shù)法演算系數(shù)表示為

(5)

分段馬斯京根法、MC演算法和MCT可變參數(shù)法的出流公式為

Ot+Δt=C1It+Δt+C2It+C3Ot

(6)

式中：Ot、Ot+Δt分別為t和t+Δt時刻下斷面出流，m3/s；It、It+Δt分別為t和t+Δt時刻的上斷面入流，m3/s。

2 研 究 流 域

沅陵站—王家河站河段地處湖南省西北部，懷化市北端，屬于沅水流域，為亞熱帶季風(fēng)氣候區(qū)，位于110°05′31″E～111°06′27″E、28°04′48″N～29°02′26″N之間，區(qū)間面積約為4 100 km2。在長期的水流侵蝕、切割和新構(gòu)造運(yùn)動的作用下，形成南北高起、中間陷落的沅陵谷地地貌。境內(nèi)山巒重疊，溪河縱橫，地形復(fù)雜。春秋兩季氣旋活動頻繁，冷暖變化大。春季及初夏多鋒面雨，夏秋之際多臺風(fēng)。沅陵縣屬于中亞熱帶季風(fēng)濕潤氣候，陽光充足，雨量充沛，四季分明，年平均氣溫16.7℃。一般8月最高，1月最低。年均降水量1 400 mm左右，常有局部暴雨出現(xiàn)，夏季降水量占全年雨量的40%。研究區(qū)域面積為6 456 km2，河段總長112 km，底坡較緩、較均勻，比降約為0.000 277，河槽糙率n為0.022～0.030。圖1為沅陵站—王家河站河段所處流域示意圖，河道徑流主要依靠上游來水，旁側(cè)只有小支流匯入，入流量均不大，因此沅陵站—王家河站河段干流受支流來水影響較小，本文不予考慮，可以看成單一河道。

圖1 研究河段示意圖Fig.1 Map of study channel

圖2 沅陵站—王家河站河段河道斷面概化Fig.2 Generalized figure of Yuanling-Wangjiahe channel section

3 參 數(shù) 率 定

3.1 河道斷面概化

根據(jù)沅陵站—王家河站河段的實(shí)測大斷面資料，考慮對河道橫斷面進(jìn)行三角形概化或梯形概化。

a.三角形概化。概化過程中保留河道水深，但會損失河道邊坡和過水面積的精度，而三角形斷面河底寬度為0 m，不符合實(shí)際斷面情況。

b.梯形概化。概化后基本與河道原斷面邊坡和過水面積一致，損失部分水深。

因此，根據(jù)河段實(shí)際斷面特性，選擇將河道斷面形狀概化為梯形。其中沅陵站河道橫斷面概化后的河底寬度為382.4 m，邊坡坡度約為20.3°。王家河站河道橫斷面概化后的河底寬度為390.5 m，邊坡坡度約為17.3°。河道概化圖見圖2。沅陵站和王家河站均在流域干流上，實(shí)際斷面情況比較接近，采用兩個站斷面概化后的算術(shù)均值作為匯流演算方法的斷面參數(shù)值，即河道河底寬度為386.4 m，邊坡坡度為18.8°。圖2中河段左岸為零點(diǎn)。

3.2 匯流演算方法參數(shù)率定

在沅陵站—王家河站河段流量資料中，選取1979—1985年共15場洪水進(jìn)行河道流量演算。其中9場用于參數(shù)率定，6場用于方法驗證。

將3種匯流演算方法的時間步長定為1 h，通過經(jīng)驗推求法，確定分段馬斯京根法參數(shù)為：河段數(shù)為10，蓄量參數(shù)k=1 h，流量比重系數(shù)ε=0.45。

MC演算法和MCT可變參數(shù)法考慮河道斷面的形狀和水力特性，在實(shí)際率定過程中，當(dāng)n分別為0.03和0.02時，后者河道匯流方法的模擬精度明顯更高。河道匯流方法的模擬精度表現(xiàn)出與河道水力因素有關(guān)，分析其原因為所選洪水中多數(shù)場次流量量級超過1萬m3/s，流量偏高造成對河道的沖刷程度更高，使得實(shí)際糙率變低，說明匯流方法考慮河道特性具有可行性。匯流方法使用牛頓法對參比流量和水位進(jìn)行迭代計算[25]，具體河道參數(shù)為：河段數(shù)15，S0=0.277×10-3，n=0.021，Δx=7 467 m，Δt=1 h。

4 結(jié)果與分析

4.1 匯流演算方法參數(shù)分析

根據(jù)沅陵站實(shí)測流量數(shù)據(jù)，利用MC演算法和MCT可變參數(shù)法計算，模擬得到王家河站的出流數(shù)據(jù)。選取模擬的末時段，分析兩種方法中的庫朗數(shù)和雷諾數(shù)在15個河段上的變化過程，以19830710號洪水過程為例，其中參數(shù)沿河段的變化曲線見圖3。

圖3 兩種匯流演算方法參數(shù)變化過程Fig.3 Parameter change process of two channel flow calculation methods for the confluence calculation

分段馬斯京根法的參數(shù)k和ε在每個分段河道中不產(chǎn)生變化，屬于線性參數(shù)。MC演算法和MCT可變參數(shù)法的演算參數(shù)與每個分段河道的參比流量Q有關(guān)。根據(jù)圖3可知，兩種方法的C和D沿河段呈曲線變化，屬于非線性參數(shù)。

4.2 河道蓄量變化結(jié)果

對于一場洪水的模擬，通常更關(guān)注模擬過程線的洪峰、洪量、峰現(xiàn)時間與實(shí)際洪水過程的誤差關(guān)系，實(shí)際上洪水的蓄量守恒狀態(tài)也對模擬效果起著一定的作用。以式(7)計算河段的蓄量變化過程。若模擬的過程線蓄量在漲落后不守恒，則必然對模擬精度存在一定程度的偏差。在保證分段馬斯京根法、MC演算法和MCT可變參數(shù)法3種方法模擬的起漲流量和完全退洪時的流量一致后，對率定期和驗證期的洪水模擬結(jié)果計算蓄量守恒誤差：

(7)

15場洪水中選取4場洪水計算匯流過程中河段蓄量的變化過程值。圖4和表1分別為分段馬斯京根法、MC演算法和MCT可變參數(shù)法模擬后河段蓄量的變化過程和蓄量相對誤差。實(shí)際模擬過程中MC演算法的蓄量變化值證明了MC演算法的蓄量不守恒。觀察3種方法在圖4的始末位置情況，可以發(fā)現(xiàn)MC演算法相較于另外兩種方法，在洪水結(jié)束時的河道蓄量值明顯高于初始值。由表1可以得到，MC演算法的蓄量相對誤差較大，最大誤差值僅略低于45%；MCT可變參數(shù)法的蓄量誤差很小，分析其原因可能為在實(shí)際計算過程中分段參數(shù)的舍入誤差。分段馬斯京根法由于匯流參數(shù)為恒定值，所以始末蓄量始終一致。

圖4 蓄量變化過程Fig.4 Change process diagram of storage volume

表1 蓄量相對誤差

4.3 總體模擬結(jié)果分析

分段馬斯京根法、MC演算法和MCT可變參數(shù)法3種方法河道匯流模擬結(jié)果見表2。選取洪峰相對誤差、洪量相對誤差、峰現(xiàn)時間誤差和確定性系數(shù)4個指標(biāo)對模擬效果進(jìn)行評價。圖5為3種方法在沅陵站—王家河站河段率定期和驗證期的總體模擬精度情況(圖中M、MC、MCT分別代表分段馬斯京根法、MC演算法和MCT可變參數(shù)法)。圖6為沅陵站—王家河站河段19820915號、19830510號、19830710號和19840607號洪水用3種方法得到的模擬流量過程線。

表2 3種方法河道匯流模擬結(jié)果

圖5 沅陵站—王家河站河段洪水模擬精度箱線圖Fig.5 Box plot of flood simulation accuracy in Yuanling-Wangjiahe channel section

圖6 沅陵站—王家河站洪水模擬過程線Fig.6 Observed and simulated flood hydrographs in Yuanling-Wangjiahe reach

4.3.1 率定期結(jié)果分析

率定期共有9場洪水，從表2可以得到分段馬斯京根法、MC演算法和MCT可變參數(shù)法的洪峰和洪量合格率均為100%，分段馬斯京根法的峰現(xiàn)時間合格率為77.8%，MC演算法和MCT可變參數(shù)法的峰現(xiàn)時間合格率均為88.9%。分段馬斯京根法、MC演算法和MCT可變參數(shù)法的確定性系數(shù)均值分別為0.870、0.889、0.891。

洪量誤差方面，19790708號洪水、19800502號洪水、19820915號洪水在3種方法的模擬下洪量誤差均超過-10%，其他場次洪量誤差均在-10%以內(nèi)；洪峰誤差方面，MC演算法和MCT可變參數(shù)法有兩場洪水洪峰誤差超過-10%，分段馬斯京根法只有一場超過-10%，分段馬斯京根法表現(xiàn)更優(yōu)；峰現(xiàn)時間誤差方面，分段馬斯京根法模擬19800518號洪水峰現(xiàn)時間誤差較大，檢查發(fā)現(xiàn)該洪水過程為多峰洪水，且兩個峰值之間相隔時間較短，分段馬斯京根法未模擬出實(shí)際最高的峰值，使得模擬峰現(xiàn)時間誤差較大。除了該次洪水，只有一場洪水峰現(xiàn)時間誤差超過10 h，3種方法很好地模擬了實(shí)際下游斷面的峰現(xiàn)時間，說明分段馬斯京根法、MC演算法、MCT可變參數(shù)法都很好地模擬出洪水在河道的擴(kuò)散和傳遞。

4.3.2 驗證期結(jié)果分析

驗證期共有6場洪水，分段馬斯京根法、MC演算法、MCT可變參數(shù)法的模擬精度結(jié)果見表2。3種方法的洪量合格率均為100%，只有19840606號洪水洪量誤差超過-10%。分段馬斯京根法的洪峰合格率為83.3%，MC演算法和MCT可變參數(shù)法的洪峰合格率略低，均為66.7%。19840606號洪水和19850703號洪水洪峰誤差相對較大，其中19840606號洪水為多峰洪水過程，3種方法在第1個洪峰處模擬較好，但在第2個洪峰處，王家河站出流量遠(yuǎn)大于沅陵站入流量，分析其原因可能為該時段流域內(nèi)降水量較大，區(qū)間入流增大。峰現(xiàn)時間誤差均控制在10 h以內(nèi)，模擬結(jié)果較好。除去19840606號洪水外，分段馬斯京根法、MC演算法和MCT可變參數(shù)法的確定性系數(shù)均超過了0.9。

率定期的9場洪水對匯流參數(shù)進(jìn)行了很好的率定，使得率定期和驗證期次洪模擬達(dá)到了較好的精度。從圖6可以看到分段馬斯京根法、MC演算法、MCT可變參數(shù)法的模擬過程線與實(shí)測比較接近，對于多峰洪水的模擬效果也比較好。

總體上看，分段馬斯京根法、MC演算法、MCT可變參數(shù)法模擬結(jié)果均比較理想。在模擬過程中，沒有考慮旁側(cè)入流，實(shí)際結(jié)果也反映出這個假設(shè)基本正確。由于MC演算法蓄量不守恒的問題，其最直接的影響是洪量模擬精度。從表2可以發(fā)現(xiàn)，除了19840430號洪水MC演算法的洪量模擬誤差低于分段馬斯京根法外，其余場次洪水洪量相對誤差均最高，印證了蓄量不守恒會增加洪量誤差的結(jié)論。從圖5可見，MC演算法和MCT可變參數(shù)法的確定性系數(shù)相對分段馬斯京根法精度更高且更為穩(wěn)定。

5 結(jié) 論

a.MC演算法和MCT可變參數(shù)法基于運(yùn)動波方程差分解的數(shù)值擴(kuò)散在一定條件下模擬擴(kuò)散波的物理擴(kuò)散，能夠反映天然河道中洪水波既傳播又?jǐn)U散的運(yùn)動規(guī)律，在實(shí)際模擬過程中考慮河道斷面形狀和水力特性，貼近實(shí)際河道，對流域內(nèi)寬淺河道的大中型洪水模擬效果較好。

b.MCT可變參數(shù)法在MC演算法的基礎(chǔ)上加入了校正系數(shù)，增強(qiáng)了參數(shù)的時變性，解決了蓄量和穩(wěn)態(tài)不一致的問題，在有限的范圍內(nèi)進(jìn)一步提高了MC演算法的穩(wěn)定性和模擬精度。

c.在河道斷面相差不大的河段中，以斷面參數(shù)平均值作為整個河段的概化形狀，結(jié)果表明模擬精度較為滿意，具有一定的可行性。

d.總體上看，分段馬斯京根法、MC演算法、MCT可變參數(shù)法在沅陵站—王家河站河段的模擬結(jié)果較好。多數(shù)洪水均為多峰洪峰，次洪中存在多個流量漲落過程，且漲落較快，3種方法對于該類洪水也具有較強(qiáng)的適應(yīng)性和可靠性。分段馬斯京根法簡便的計算步驟仍然可以帶來良好的模擬效果。在實(shí)際模擬中，MC演算法和MCT可變參數(shù)法與分段馬斯京根法取得的效果不相上下，甚至更優(yōu)。