李致家,何 蒙,閆鳳翔,胡友兵,劉志雨,童冰星

(1.河海大學(xué)水文水資源學(xué)院,江蘇 南京 210098; 2.河北省水文水資源勘測(cè)局,河北 石家莊 050031;3.淮河水利委員會(huì)水文局,安徽 蚌埠 233000; 4.水利部水文局,北京 100053)

根據(jù)河段上斷面入流過(guò)程,應(yīng)用洪水波運(yùn)動(dòng)的數(shù)學(xué)描述方法解出流量的時(shí)空變化函數(shù)或者某一指定下斷面的流量過(guò)程,稱(chēng)為河道洪水演算[1]。常用的演算方法有2 種：基于水量平衡與槽蓄方程的水文學(xué)方法和基于圣維南方程組及其簡(jiǎn)化算法的水力學(xué)方法[2]。水文學(xué)方法計(jì)算相對(duì)簡(jiǎn)單、資料需求相對(duì)較低，易于在實(shí)時(shí)洪水預(yù)報(bào)中應(yīng)用,在我國(guó)洪水預(yù)報(bào)中應(yīng)用廣泛,其中最為典型的方法是馬斯京根法[1,3-4]。馬斯京根法假定出流與槽蓄量成線性函數(shù)關(guān)系,同時(shí)流量在計(jì)算時(shí)段內(nèi)和沿程變化呈直線分布,參數(shù)要利用實(shí)測(cè)洪水資料率定,或者根據(jù)參數(shù)的水力學(xué)公式估算[3]。我國(guó)河流眾多,水流條件復(fù)雜多變,同時(shí)受人類(lèi)活動(dòng)影響,河道的斷面形狀和水力特性可能發(fā)生變化，因此，本文選用Muskingum-Cunge-Todini(MCT)可變參數(shù)法和非線性水庫(kù)法進(jìn)行河道洪水演算。選用的2種演算方法能夠考慮河道斷面形狀和水力特性,在國(guó)外已經(jīng)取得較為廣泛的應(yīng)用[5-6]。選取濕潤(rùn)地區(qū)淮河中游吳家渡—小柳巷河段及半干旱地區(qū)滹沱河中游黃壁莊—北中山河段為典型單一河段驗(yàn)證這2種演算方法在不同在流域的適用性,并將模擬結(jié)果與馬斯京根法進(jìn)行分析比較。

1 方 法 原 理

1.1 馬斯京根法

19世紀(jì)70年代圣維南建立明渠非恒定流偏微分方程組，為洪水研究奠定了理論基礎(chǔ)[7]。McCarthy于1938年提出馬斯京根法，開(kāi)創(chuàng)了洪水演算的水文學(xué)方法[8]，其流量演算方程式為

O2=C0I2+C1I1+C2O1

(1)

式中：O1、O2——計(jì)算時(shí)段始、末的河段出流量；I1、I2——計(jì)算時(shí)段始、末的河段入流量；C0、C1、C2——流量演算系數(shù)；K——蓄量常數(shù)；x——流量比重系數(shù)，與洪水波的坦化變形程度有關(guān)；Δt——計(jì)算時(shí)段長(zhǎng)。

馬斯京根法采用線性有限差解，要求水面線為直線且流量沿程線性變化，因此演算時(shí)段Δt應(yīng)等于或接近K值[9]。但在長(zhǎng)河段的情況下，K遠(yuǎn)大于Δt，這時(shí)需要采用分段馬斯京根法[10]。

1.2 MCT可變參數(shù)法

1969年Cunge在Muskingum法的基礎(chǔ)上提出M-C方法，開(kāi)辟了應(yīng)用運(yùn)動(dòng)波數(shù)值擴(kuò)散進(jìn)行洪水演算的新途徑[11]。1978年P(guān)ouce和Yevjevich在Cunge的基礎(chǔ)上，考慮到波速的變動(dòng)問(wèn)題，建議采用變動(dòng)參數(shù)表示動(dòng)力波速和水面寬度，提出了變動(dòng)參數(shù)M-C方法[12]。但在變動(dòng)參數(shù)M-C方法的實(shí)際應(yīng)用中，發(fā)現(xiàn)了2個(gè)問(wèn)題：一是質(zhì)量不守恒問(wèn)題[13]；二是蓄量不穩(wěn)定問(wèn)題[14]。由此，Todini對(duì)變動(dòng)參數(shù)M-C方法進(jìn)行修正，提出一種新的質(zhì)量守恒和穩(wěn)定狀態(tài)協(xié)調(diào)的Muskingum-Cunge-Todini可變參數(shù)法(簡(jiǎn)稱(chēng)為MCT可變參數(shù)法)(式(2))，并用于TOPKAPI模型中[14-17]。

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

(2)

1.3 非線性水庫(kù)法

運(yùn)動(dòng)學(xué)方法中采用曼寧公式近似動(dòng)量方程，得到非線性水庫(kù)方程的一般形式：

(3)

式中：y——河道蓄水量；a、b、c——系數(shù)，在每個(gè)時(shí)間步長(zhǎng)內(nèi)為不變量。

劉志雨等[18-19]、Todini等[20]發(fā)現(xiàn)對(duì)式(3)中的指數(shù)項(xiàng)進(jìn)行二次方程近似，可將復(fù)雜的指數(shù)方程轉(zhuǎn)化為低階方程并得到解析解。假設(shè)河道為樹(shù)形結(jié)構(gòu)，河段為三角形斷面，非線性水庫(kù)方程的具體表達(dá)式為

(4)

式中：Vc——河道蓄水量；rc——旁側(cè)入流量，包括地表徑流量和到達(dá)河道的土壤排水量；Qc——來(lái)自河道上游的入流量；S0——河底比降，假定與地形坡度相同；n——河道曼寧糙率系數(shù)；γ——河道側(cè)邊與水平線的夾角；Δx——計(jì)算步長(zhǎng)。

表2 2005年、2007年吳家渡—小柳巷河段慣性項(xiàng)

注:Sl為局地慣性項(xiàng),Sx為遷移慣性項(xiàng)。

2 淮 河 算 例

2.1 河道概況

淮河干流中游吳家渡—小柳巷河段長(zhǎng)約94 km,河床糙率約為0.02～0.025[21]。從吳家渡到小柳巷依次分布有方邱湖行洪區(qū)、臨北縷堤湖行洪區(qū)、花園湖行洪區(qū)、香浮段行洪區(qū)、潘村段行洪區(qū)[22]。除花園湖行洪區(qū)偶爾使用之外,受上游諸多閘門(mén)調(diào)度與行蓄洪區(qū)攔蓄作用影響,其他行洪區(qū)很少有機(jī)會(huì)使用。由于北淝河上的沫河口閘、懷洪新河上的五河閘常年關(guān)閉,這2條河的來(lái)水少有機(jī)會(huì)匯入干流,因此吳家渡—小柳巷河段干流受支流來(lái)水影響較小,本文不予考慮。另外由于小柳巷以下河道較為順直,且離洪澤湖較遠(yuǎn),因此受洪水頂托作用較小。

吳家渡到小柳巷縱剖面變化劇烈,最淺處與最深處深泓高程相差較大[22]?？v剖面的急劇變化導(dǎo)致水流更加紊亂,河道形狀阻力加大。根據(jù)吳家渡、臨淮關(guān)、五河、小柳巷4個(gè)站的實(shí)測(cè)斷面資料,計(jì)算出平均河底比降約為0.002 06%,并將河道斷面形狀概化為河床夾角沿程變化的三角形,體現(xiàn)河段上下游的斷面形狀變化。

2.2 洪水波水力特征分析

2.2.1 附加比降分析

直接測(cè)定附加比降十分困難,所以一般認(rèn)為附加比降SΔ近似等于水面比降S與河底比降S0之差[3]。選取2003年、2005年、2007年作為典型年,根據(jù)實(shí)測(cè)的吳家渡、臨淮關(guān)、五河、小柳巷的日均水位資料, 選擇吳家渡、小柳巷兩站水勢(shì)平穩(wěn)及流量接近的日期計(jì)算吳家渡—小柳巷河段水面比降S,并推求出洪水波的附加比降,結(jié)果見(jiàn)表1。

表1 2003年、2005年、2007年吳家渡—小柳巷河段水面比降和附加比降

注:各河段的河底比降S0均為0.002 06%,SΔ為負(fù)值是因?yàn)橛?jì)算日期處于落洪期。

一般約定波前即漲洪段的附加比降為正,波后即落洪段的附加比降為負(fù)。2003-09-27、2005-10-09臨淮關(guān)—五河段的附加比降為正值,是因?yàn)?003年、2005年洪水較大,臨淮關(guān)—五河段的水面比降大于吳家渡—臨淮關(guān)段和五河—小柳巷段,是阻水河段,使得洪水波的附加比降出現(xiàn)正值?？梢钥闯?落洪期附加比降占河底比降為4%～10%,且漲洪期附加比降一般大于落洪期附加比降的絕對(duì)值,占比將更大,其影響不可忽視。

2.2.2 慣性項(xiàng)分析

由于2003年沒(méi)有實(shí)測(cè)流速資料,僅對(duì)吳家渡—小柳巷河段2005年、2007年計(jì)算時(shí)段內(nèi)的慣性項(xiàng)進(jìn)行分析,結(jié)果見(jiàn)表2。本河段洪水波慣性項(xiàng)僅占河底比降1%～4%,故可以忽略。

由表1～2可知,動(dòng)力方程中慣性項(xiàng)占河底比降的比例相當(dāng)小,可以忽略;而附加比降項(xiàng)占河底比降的比例較大,不能忽略。因此吳家渡—小柳巷河段擴(kuò)散波特征明顯。

2.3 河道洪水演算方法的參數(shù)率定

選取2003—2017年大、中、小洪水共15場(chǎng)進(jìn)行洪水流量過(guò)程模擬,前9場(chǎng)用于參數(shù)率定,其余6場(chǎng)用于檢驗(yàn)?zāi)Ｐ?。為減少土壤、地表和河道等蓄水初始條件對(duì)計(jì)算精度的影響,將洪水發(fā)生的前一個(gè)月作為預(yù)熱期,以吳家渡實(shí)測(cè)流量作為輸入進(jìn)行河道洪水演算。通過(guò)人工試錯(cuò),結(jié)合流域?qū)嶋H情況,將3種演算方法的時(shí)間步長(zhǎng)定為1 h,河段分為16段,確定河床糙率為0.022,分段馬斯京根法x=-0.3。

2.4 河道洪水演算結(jié)果

由表3可知:(a)非線性水庫(kù)法和分段馬斯京根法的洪量相對(duì)誤差都在±10%以?xún)?nèi),MCT可變參數(shù)法洪量相對(duì)誤差只有1場(chǎng)超過(guò)-10%,其余場(chǎng)次結(jié)果均較好。(b)非線性水庫(kù)法和分段馬斯京根法洪峰相對(duì)誤差都有2場(chǎng)超過(guò)±20%,而MCT可變參數(shù)法的洪峰相對(duì)誤差都在±20%以?xún)?nèi),精度高于其他2種方法。(c)3種演算方法的峰現(xiàn)時(shí)間誤差在各場(chǎng)次洪水間變化較大,但是平均誤差都在±5 h以?xún)?nèi)?？傮w上來(lái)看MCT可變參數(shù)法的峰現(xiàn)時(shí)間誤差除在20050707、20160616、20170923場(chǎng)次洪水模擬中較大,其余場(chǎng)次模擬結(jié)果較好。(d)非線性水庫(kù)法和MCT可變參數(shù)法的確定性系數(shù)有2場(chǎng)小于0.9,分段馬斯京根法的確定性系數(shù)有1場(chǎng)小于0.9。

表3 淮河3種河道洪水演算方法模擬結(jié)果

注:EW為洪量相對(duì)誤差,EQ為洪峰相對(duì)誤差,ET為峰現(xiàn)時(shí)間誤差,DC為確定性系數(shù)。

圖1 小柳巷洪水模擬結(jié)果Fig.1 Flood simulation results in Xiaoliuxiang

綜合表3和圖1分析,可得:MCT可變參數(shù)法對(duì)20130825場(chǎng)次洪水模擬結(jié)果較差,洪量相對(duì)誤差為-11%,洪峰相對(duì)誤差為-16.54%。是因?yàn)楸緢?chǎng)洪水量級(jí)明顯小于率定期其他場(chǎng)次洪水,區(qū)間入流影響較大且洪水歷時(shí)偏短,洪水陡漲陡落。由于偽數(shù)值擴(kuò)散的存在,MCT可變參數(shù)法在小洪水模擬過(guò)程中擴(kuò)散效果偏大,導(dǎo)致洪峰流量偏小,洪量偏小。而基于水量平衡方程的非線性水庫(kù)法和分段馬斯京根法模擬效果較好。個(gè)別場(chǎng)次洪水MCT可變參數(shù)法的峰現(xiàn)時(shí)間模擬結(jié)果較差的原因是:(a)本次參數(shù)是人工率定結(jié)果,率定過(guò)程中側(cè)重于確定性系數(shù)和洪峰相對(duì)誤差,不能全面考慮洪水過(guò)程的特征指標(biāo)。(b)河道斷面概化較簡(jiǎn)單,沒(méi)有考慮灘地的行洪過(guò)程。(c)河床糙率對(duì)峰現(xiàn)時(shí)間的影響。參數(shù)率定中河床糙率對(duì)洪峰流量、峰現(xiàn)時(shí)間都有影響,改變峰現(xiàn)時(shí)間的同時(shí)可能會(huì)影響模擬洪峰值,使得峰現(xiàn)時(shí)間的模擬結(jié)果不理想。(d)淮河中游河段倒比降的影響。根據(jù)實(shí)測(cè)資料及相關(guān)研究[23-24],淮河中游以浮山為轉(zhuǎn)折點(diǎn),吳家渡—浮山段平均深泓沿程降低,呈正比降;浮山—小柳巷段平均深泓沿程增大,呈負(fù)比降,致使淮河中游洪水無(wú)法通暢下泄,模擬結(jié)果與實(shí)際情況相差較大。

總體上看3種演算方法在研究河段上都取得了較好的模擬效果,且MCT可變參數(shù)法和非線性水庫(kù)法在研究河道上模擬精度較高,洪峰合格率均在86%以上,確定性系數(shù)均大于0.8,表明這2種方法可以為河道洪水演算提供新的解決方案。

3 滹沱河算例

3.1 河道概況

滹沱河為海河流域子牙河系兩大支流之一,屬于季節(jié)性河流,河道在枯季常有干涸斷流現(xiàn)象,僅在暴雨期間形成河道徑流。河床以中細(xì)砂為主,滲透性強(qiáng)。黃壁莊水庫(kù)—北中山河段為滹沱河的主干流,河長(zhǎng)約110 km,平均河底比降約為0.051 6%。流經(jīng)平原,河道寬淺,可以概化為矩形河道。本河段徑流主要靠黃壁莊水庫(kù)泄洪,區(qū)間入流可忽略,以河道洪水演算為主。

傳統(tǒng)的馬斯京根河道流量演算法基于水量平衡對(duì)河道洪水進(jìn)行演算,但對(duì)于常年斷流、河水與地下水長(zhǎng)期處于脫節(jié)狀態(tài)的河道,洪水傳播過(guò)程中下滲劇烈,滲漏損失量很大,天然狀態(tài)下洪水演進(jìn)過(guò)程發(fā)生顯著變化,與模型基本原理不相符[25-26]。所以,考慮滲漏損失的洪水演進(jìn)模擬成為半干旱地區(qū)洪水模擬的關(guān)鍵問(wèn)題。

3.2 考慮河道滲漏的河道演算方法

目前,河道滲漏損失估算及模擬方面研究已取得了一定成果,國(guó)內(nèi)研究多是直接采用下滲公式方法[27-30]。本文借鑒考慮河道滲漏的分段馬斯京根法,對(duì)MCT可變參數(shù)法和非線性水庫(kù)法增加河道滲漏計(jì)算。采用霍頓下滲公式(式(5))計(jì)算每段的下滲率,轉(zhuǎn)化成下滲流量。在入流中扣除后得到每段的凈入流量,然后演算到流域出口斷面。

ft=(f0-fc)ekt+fc

(5)

式中:ft——下滲率;fc——穩(wěn)定下滲率;f0——初始下滲率;k——與擴(kuò)散率有關(guān)的系數(shù),與河床的物理特性有關(guān);t——下滲歷時(shí)。

3.3 河道洪水演算方法的參數(shù)率定

由于實(shí)測(cè)洪水資料比較少,選取研究河段20世紀(jì)80年代以后2場(chǎng)典型洪水19880805、19960831進(jìn)行模擬,探討考慮河道滲漏的MCT可變參數(shù)法和非線性水庫(kù)法在半干旱區(qū)的適用性。洪水的傳播時(shí)間、輸水損失與流量大小、河床干旱程度有較密切的關(guān)系。首先根據(jù)河道前期干濕狀況、河道特性等河段具體情況確定一組模型參數(shù),得出河道下游出流過(guò)程,將此過(guò)程與河道下游實(shí)際出流過(guò)程相比較,同時(shí)參考下游洪峰流量及起漲點(diǎn)的擬合情況進(jìn)行綜合分析。經(jīng)過(guò)率定,下滲參數(shù)fc取5.2 mm/h,k取0.2,f0取1 500 mm/h。將3種演算方法的時(shí)間步長(zhǎng)定為1 h,河段分為16段,確定河床糙率為0.03,分段馬斯京根法x取0.029。

3.4 河道洪水演算結(jié)果分析

從表4可以看出:(a)19880805場(chǎng)次洪水3種演算方法的洪峰相對(duì)誤差均在±20%以?xún)?nèi),但都明顯偏小;峰現(xiàn)時(shí)間模擬結(jié)果均較好;確定性系數(shù)整體不高。(b)19960831場(chǎng)次洪水的3種演算方法模擬結(jié)果均良好。3種方法的洪峰相對(duì)誤差均在±2%左右;峰現(xiàn)時(shí)間模擬結(jié)果均較好;確定性系數(shù)都大于0.9;整體上MCT可變參數(shù)法的精度高于其他2種方法。圖2為北中山洪水模擬結(jié)果。

表4 滹沱河3種河道洪水演算方法模擬結(jié)果

圖2 北中山洪水模擬結(jié)果Fig.2 Flood simulation results in Beizhongshan

19880805場(chǎng)次洪水模擬結(jié)果不好的原因可能有:(a)本場(chǎng)洪水是在河道連續(xù)8 a無(wú)水的情況下,黃壁莊水庫(kù)上發(fā)生洪水開(kāi)始放水,河道見(jiàn)水歷時(shí)長(zhǎng)達(dá)180 h左右,河道滲漏損失水量約占來(lái)水量的52%,輸水條件較為復(fù)雜,下滲過(guò)程模擬比較困難。(b)參數(shù)率定過(guò)程中側(cè)重于洪水起漲部分的模擬效果,洪峰及落水段沒(méi)有兼顧好。(c)洪水量級(jí)較小,與19960831場(chǎng)次洪水差別較大,且2場(chǎng)洪水時(shí)間間隔較長(zhǎng),下墊面因素可能發(fā)生變化,下滲條件有所差別,使用同一套參數(shù)模擬結(jié)果較差。

4 結(jié) 語(yǔ)

a. MCT可變參數(shù)法基于運(yùn)動(dòng)波方程差分解的數(shù)值擴(kuò)散在一定條件下模擬擴(kuò)散波的物理擴(kuò)散,能夠反映天然河道中洪水波既傳播又?jǐn)U散的運(yùn)動(dòng)規(guī)律,并引入修正因子保證計(jì)算過(guò)程中質(zhì)量守恒和蓄量穩(wěn)定。本方法考慮了河道斷面形狀和水力學(xué)糙率,參數(shù)具有明確的物理意義,對(duì)平原緩坡河段的大中型洪水也適用。

b. 非線性水庫(kù)法在計(jì)算過(guò)程中考慮了河道斷面情況和水力學(xué)糙率,在研究河段模擬效果較好,精度和分段馬斯京根法相近,具有較強(qiáng)的適用性。

c. 分段馬斯京根法計(jì)算簡(jiǎn)便,經(jīng)驗(yàn)性較強(qiáng),參數(shù)具有一定的物理意義,在研究河段模擬效果較好,具有廣泛適用性。

整體來(lái)看3種河道洪水演算方法在吳家渡—小柳巷河段的模擬結(jié)果均比較滿(mǎn)意,精度差異不大。但是,MCT可變參數(shù)法和非線性水庫(kù)法的參數(shù),可以根據(jù)河道斷面的大斷面資料和水力學(xué)特性直接估算,考慮河道斷面變化對(duì)洪水的影響,可以為河道洪水預(yù)報(bào)提供更多的選擇,拓廣了現(xiàn)有河道洪水預(yù)報(bào)方法,同時(shí)也為相對(duì)復(fù)雜的河道洪水預(yù)報(bào)問(wèn)題研究提供了一定的參考借鑒。同時(shí)本文也在北方半干旱地區(qū)初步分析了考慮河道滲漏的3種河道演算方法,模擬結(jié)果表明研究思路和方法具有一定的合理性和可行性。鑒于河道洪水預(yù)報(bào)的復(fù)雜性,方法對(duì)其他流域的適用性還需進(jìn)一步論證和檢驗(yàn)。