王愷禎，王 軍，隋覺義

（1.合肥工業(yè)大學(xué) 土木與水利工程學(xué)院，安徽 合肥 230009；2.北英屬哥倫比亞大學(xué) 環(huán)境工程系，加拿大 喬治王子城 V2N 4Z9）

1 研究背景

特殊的地理位置和氣候條件使黃河寧蒙河段每年凌汛期（當(dāng)年11月至次年3月）發(fā)生不同程度的凌情［1］，冰凌會危及橋梁、涵洞、河道工程和其他水工建筑物的安全［2-3］。在封、開河時(shí)期，如果冰凌嚴(yán)重堵塞河道，造成水位大幅度地抬高而引發(fā)漫灘或決堤，會給黃河灘區(qū)和兩岸人民群眾的生命財(cái)產(chǎn)安全帶來威脅和損失，準(zhǔn)確及時(shí)的凌汛預(yù)報(bào)能夠?yàn)榉姥垂ぷ魈峁Q策支持。

處于穩(wěn)定流狀態(tài)的河道水流，由于暴雨徑流、水電站運(yùn)行或閘壩放水等原因，使原來的水面受到干擾而形成不穩(wěn)定波動，即洪水波，波流量從河道上游向下游運(yùn)動，就產(chǎn)生了河道里的洪水波運(yùn)動。洪水波在演進(jìn)過程中會發(fā)生變形［4］，生產(chǎn)實(shí)際中，洪水波的傳播機(jī)制是洪水預(yù)報(bào)、防洪規(guī)劃和水利工程防洪控制運(yùn)用中重要的組成部分，其研究具有重要的理論和實(shí)際意義。

目前廣泛應(yīng)用于洪水波的傳播計(jì)算的馬斯京根法由McCarthy于1938年應(yīng)用于美國馬斯京根河流域而得名，后經(jīng)一系列研究，到1978年P(guān)ouce和Yevjevich將變動水面寬度及動力波速度的變動參數(shù)應(yīng)用于1969年Cunge提出的馬斯京根-康吉模型［5］上，馬斯京根模型此時(shí)才基本完善。該方法在實(shí)際應(yīng)用中的一個(gè)重要問題是模型參數(shù)的估計(jì)。在傳統(tǒng)方法中，利用試算法對參數(shù)K、x進(jìn)行率定，再計(jì)算出決策變量C0、C1、C2，但是試算法具有盲目性、不確定性等缺陷，因此有學(xué)者選擇了直接優(yōu)化決策變量C0、C1、C2，得到較好的效果，主要方法有最小二乘法［6］，遺傳算法［7］，粒子群算法［8］等。馬斯京根法的參數(shù)K反映洪水波的傳播時(shí)間，而x則反映河段蓄量的大小及調(diào)蓄能力。

對于洪水演進(jìn)，可采用水力學(xué)的非恒定流法和水文學(xué)的馬斯京根法來計(jì)算。當(dāng)河流有冰存在時(shí)，由于冰水存在能量、質(zhì)量交換，使得計(jì)算難度加大，而馬斯京根法可以從水量平衡方程中考慮到冰水量的相互轉(zhuǎn)化問題，本文擬應(yīng)用馬斯京根法對冰期冰蓋下洪水波的變形進(jìn)行嘗試性研究計(jì)算。

2 考慮有冰存在的馬斯京根法

在水文循環(huán)過程中，任一區(qū)域、任一時(shí)段內(nèi)輸入水量與輸出水量之差等于其蓄水量的變化量，根據(jù)質(zhì)量守恒定律可寫出某一區(qū)域任一時(shí)段（Δt）內(nèi)的水量平衡方程：

冰期進(jìn)行水量平衡計(jì)算時(shí)，若忽略上游水流向下游傳播過程中轉(zhuǎn)化為冰的水量，以及上游冰蓋或浮冰消融流入下游水流中的水量，而完全使用水量平衡方程的一般公式，則會增大下游斷面計(jì)算流量的誤差，考慮冰量存在的水量平衡方程應(yīng)為：

式中：I1、I2為上游斷面時(shí)段始、末流量，m3/s；Q1、Q2為下游斷面時(shí)段始、末流量，m3/s；Δt為所取的計(jì)算時(shí)長，s；ΔW為時(shí)段槽蓄水量的變化量，m3；Δd1、Δd2為上、下游斷面冰蓋厚度的變化值，m；ΔB1、ΔB2為上、下游斷面冰蓋寬度的變化值，m；vˉ冰1、vˉ冰2為上、下游斷面冰蓋下的平均流速，m/s。

天然河道洪水波運(yùn)動的動力方程為：

式中：Z為水位，m；L為河段長，m；v為斷面平均流速，m/s；g為重力加速度，m/s2；為水面比降；Sf為摩阻項(xiàng)，用曼寧公式計(jì)算；為慣性項(xiàng)。式（3）是以牛頓第二定律為基礎(chǔ)建立起來的，反應(yīng)了能量守恒，對于一般的天然河道水流，動力方程中慣性項(xiàng)較其他項(xiàng)要小兩個(gè)數(shù)量級，通常忽略。對無旁側(cè)入流的冰期河段，忽略慣性項(xiàng)，式（3）可寫成：

式中：vc為冰蓋下過水?dāng)嗝娴钠骄魉?，m/s；C為謝才系數(shù)；Rc為冰蓋下過水?dāng)嗝娴乃Π霃?，m；nc為冰蓋下的綜合糙率，其計(jì)算方法參見文獻(xiàn)［9］。

冰蓋生成時(shí)，水流的流動結(jié)構(gòu)發(fā)生顯著的改變，一般以最大流速線將流動沿水深方向分為冰蓋區(qū)和床面區(qū)兩層，對于概化的矩形斷面，床面區(qū)濕周χb=b+2H，冰蓋區(qū)濕周χi=b：

冰期河槽蓄水量Wc為

式中：H為冰蓋下流動的水深，m；b為斷面寬，m；Ac為冰蓋下過水?dāng)嗝婷娣e，m2；為冰蓋下平均過水?dāng)嗝婷娣e，m2；為冰蓋下平均流動的水深，m。

將式（8）代入到式（7）后得）

可寫成函數(shù)形式

假定冰期槽蓄水量Wc和入流量Ic、出流量Qc存在線性關(guān)系，馬斯京根法在忽略慣性項(xiàng)下的槽蓄方程為：

式中：Jc為冰期斷面的水力坡度；K和x為馬斯京根法的流量演算參數(shù)，K值等于在相應(yīng)蓄量下恒定流狀態(tài)的河段傳播時(shí)間，K值越小，洪水波傳播時(shí)間越短，反之越長；x是反應(yīng)河槽調(diào)節(jié)作用的一個(gè)指標(biāo)，即反應(yīng)洪水波傳播過程變形的程度，x減小，洪水波變形增加，反之減少。

進(jìn)一步分析K、x和糙率的關(guān)系，K的計(jì)算公式為：

明流時(shí)期，概化的矩形斷面洪水波波速c與斷面平均流速v的關(guān)系為［10］：

由曼寧公式及式（13）、式（14）可以推出明流時(shí)期K的計(jì)算公式為：

x的計(jì)算公式為［11］：

穩(wěn)定流流量Q0的計(jì)算公式為［12］

由曼寧公式及式（13）、式（14）、式（16）可以推出明流時(shí)期x的計(jì)算公式為：

式中：H0為相應(yīng)的穩(wěn)定流水位，m；J為明流期時(shí)斷面的水力坡度；n為河床糙率系數(shù)；Qb為基流流量，m3/s；Qp為洪峰流量，m3/s。

由式（15）、式（18）可知，K值與糙率成正比，x與糙率成反比，即糙率越大，洪水波傳播時(shí)間越長，洪水波的變形越大。

3 寧蒙河段冰期馬法計(jì)算驗(yàn)證

如圖1所示，寧蒙河段地處黃河流域最北端，從寧夏的石嘴山市和內(nèi)蒙古伊克昭盟的拉僧廟入境，至伊克昭盟準(zhǔn)格爾旗馬柵的榆樹灣出境，干流全長840 km，總落差162.5 m。河道流向從石嘴山至巴彥高勒為由西南流向東北，巴彥高勒至包頭市為自西向東，包頭市至清水河縣的喇嘛灣為由西北流向東南，喇嘛灣至出境為由北向南［13］。

圖1 寧蒙河段示意圖

圖2 開河期流量特性

選擇黃河寧蒙河段的巴彥高勒（東經(jīng)106°02'，北緯40°19'）、三湖河口（東經(jīng)106°47'，北緯40°36'）、包頭（東經(jīng)109°55'，北緯40°32'）、頭道拐（東經(jīng)111°04'，北緯40°16'）為研究的水文測站，采用2014—2016年流量數(shù)據(jù)進(jìn)行計(jì)算分析。

圖2（a）為2015年包頭和頭道拐水文站開河期的流量過程線，下游頭道拐水文站比上游包頭水文站提前4天半開河，其流量過程出現(xiàn)了兩場比較大的洪水，前一場洪水由下游提前開河時(shí)槽蓄水量的釋放形成，洪水呈“尖瘦”狀；后一場洪水由上游開河，槽蓄水量向下游傳播形成，洪水形狀相對前一場較“矮胖”，洪量更大。圖2（b）為2016年包頭和頭道拐水文站開河期的流量過程線，上游包頭水文站先于下游頭道拐水文站開河，洪峰演進(jìn)過程沿程增大，上下游洪峰間隔時(shí)間沿程減小。

選取穩(wěn)封期到開河期前夕、開河期、明流期三個(gè)計(jì)算時(shí)期，應(yīng)用最小二乘原理進(jìn)行馬斯京根法參數(shù)K、x的最優(yōu)估計(jì)，并對其計(jì)算值的可靠性進(jìn)行驗(yàn)證，對流量過程、洪峰流量及洪量進(jìn)行誤差計(jì)算。限于篇幅，對于流量過程的演算僅列出三湖河口-包頭河段，如圖3所示。

圖3 包頭斷面演算流量與實(shí)際流量過程對比

從數(shù)值上分析流量誤差，穩(wěn)封期到開河期前夕的誤差小于明流期的誤差小于開河期的誤差，計(jì)算值與實(shí)測值的數(shù)值和趨勢皆比較接近，流量演算擬合程度高，洪峰流量和洪量的誤差及確定性系數(shù)的計(jì)算結(jié)果如表1所示。

從確定性系數(shù)上看，馬斯京根法計(jì)算下的明流時(shí)期高于有冰存在的兩種情況，有冰存在的確定性系數(shù)最小為0.793，仍比較高，故認(rèn)為將馬斯京根法應(yīng)用于冰期是可行的。

表1 洪峰流量及洪量的誤差分析

4 封凍期洪水波傳播變形分析

將封凍期冰蓋的發(fā)展過程分為兩部分，以冰蓋達(dá)到最厚為臨界點(diǎn)，從初封到冰蓋達(dá)到最厚視為冰蓋的凍結(jié)增厚過程，冰蓋從最厚到瓦解消融視為冰蓋的融化減薄過程。參數(shù)x能夠間接反映洪水波傳播的變形，對任一時(shí)段，可通過馬斯京根方程求出參數(shù)x值，以及同時(shí)段的平均冰厚值d，從而建立參數(shù)x與平均冰厚d的散點(diǎn)關(guān)系圖。

根據(jù)2013—2017年封凍期的流量及冰厚數(shù)據(jù)進(jìn)行計(jì)算。每隔5～6 d建立一組馬斯京根法的流量方程，一組方程能夠求解得到一個(gè)x值。計(jì)算時(shí)對流量進(jìn)行內(nèi)插，時(shí)間間隔取Δt=8h，冰厚d取時(shí)段始末的均值。冰蓋凍結(jié)增厚時(shí)期參數(shù)x與平均冰厚d的散點(diǎn)關(guān)系圖如圖4所示。

圖4 冰蓋凍結(jié)增厚過程參數(shù)x與冰厚d關(guān)系

冰蓋的凍結(jié)增厚過程對洪水波變形的影響可分為3個(gè)階段。圖4（a）為第一個(gè)階段，冰蓋處于初封期，冰蓋的增厚使洪水波變形增大，這一階段冰蓋增厚速度在0.01～0.03 m/d，巴彥高勒-三湖河口段的計(jì)算冰厚范圍在0.25～0.33m，三湖河口—包頭段的計(jì)算冰厚范圍在0.14～0.27 m，包頭-頭道拐段的計(jì)算冰厚范圍在0.20～0.34 m；之后進(jìn)入第二階段，如圖4（b）所示，冰蓋的增厚使洪水波變形減小，這一階段冰蓋增厚的速度穩(wěn)定且緩慢，增厚速度小于0.01 m/d，巴彥高勒-三湖河口段的計(jì)算冰厚范圍在0.33～0.42 m，三湖河口-包頭段的計(jì)算冰厚范圍在0.27～0.43 m，包頭-頭道拐段的計(jì)算冰厚范圍在0.34～0.47m；到了第三階段，如圖4（c）所示，洪水波的變形又相對增大，這一階段冰厚幾乎沒有變化，而流量持續(xù)增大，氣溫逐漸上升，有些天數(shù)日最大氣溫已達(dá)0℃以上，巴彥高勒-三湖河口段的計(jì)算冰厚范圍在0.42～0.58 m，三湖河口-包頭段的計(jì)算冰厚范圍在0.43～0.62 m，包頭-頭道拐段的計(jì)算冰厚范圍在0.47～0.60 m。

綜上所述，在冰厚快速增加的第一階段，冰蓋的出現(xiàn)增大了過水?dāng)嗝娴臐裰芎妥枇Γ?4］，河段綜合糙率增加，洪水波變形增大；隨后在冰厚緩慢增加的第二階段，由于水力磨蝕作用，冰蓋底部漸趨光滑，糙率降低，洪水波變形減?。蛔詈笤诒駧缀醪话l(fā)生改變但流量增大的第三階段，由于冰厚保持最大值狀態(tài)，穩(wěn)封期綜合糙率沒有大的改變，由式（18）知，Q0與x成反比，流量增大，通過影響Qp，使Q0增大，x減小，故洪水波變形增大。

冰蓋融化減薄期參數(shù)x與平均冰厚d的關(guān)系圖如圖5所示。

圖5 冰蓋融化減薄過程參數(shù)x與冰厚d關(guān)系

冰蓋融化減薄的時(shí)期，由于冰蓋融化至0.35 m以下的冰厚數(shù)據(jù)極少，故計(jì)算范圍有限。計(jì)算范圍內(nèi)的洪水波變形隨冰蓋的消融而增大，巴彥高勒—三湖河口段的計(jì)算冰厚范圍為0.35～0.57 m；三湖河口—包頭段的計(jì)算冰厚范圍為0.22～0.6 m；包頭—頭道拐段的計(jì)算冰厚范圍為0.36～0.52 m。

將冰蓋凍結(jié)增厚和冰蓋融化減薄過程對比分析，巴彥高勒—三湖河口段的計(jì)算冰厚范圍在0.35～0.58 m，三湖河口—包頭段的計(jì)算冰厚范圍在0.22～0.60 m，包頭—頭道拐段的計(jì)算冰厚范圍在0.34～0.52 m。計(jì)算結(jié)果如圖6所示。

圖6 冰蓋凍結(jié)增厚與融化減薄過程參數(shù)x與冰厚d關(guān)系對比

冰蓋凍結(jié)增厚過程和冰蓋融化減薄過程下洪水波的變形程度并不一致，且兩個(gè)過程的參數(shù)x均會在冰蓋凍結(jié)增厚對洪水波傳播變形的第三階段有交點(diǎn)，交點(diǎn)之前的冰厚范圍，冰蓋融化減薄時(shí)的洪水波變形更大，交點(diǎn)之后的冰厚范圍，冰蓋凍結(jié)增厚時(shí)的洪水波變形更大。冰厚的計(jì)算范圍如表2所示。

由于洪水波傳播過程的變形程度為水深和冰蓋糙率的綜合作用，冰蓋融化減薄時(shí)的洪水波變形更多，說明該冰厚范圍內(nèi)，相同冰厚下冰蓋融化減薄過程的冰下過流能力更??；冰蓋凍結(jié)增厚時(shí)的洪水波變形更多，說明該冰厚范圍內(nèi)，相同冰厚下冰蓋融化減薄過程的冰下過流能力更大。

表2 各河段冰厚與洪水波變形分析

5 結(jié)論

基于馬斯京根法，選取黃河寧蒙河段為研究對象，經(jīng)過分析，得出以下結(jié)論：

通過K、x的計(jì)算公式，分析其與糙率的關(guān)系，得出K值與糙率成正比，x與糙率成反比，即糙率越大，洪水波傳播時(shí)間越長，洪水波的變形越大。

對計(jì)算出的馬斯京根法參數(shù)K、x的可靠性進(jìn)行驗(yàn)證，得出計(jì)算值與實(shí)測值的數(shù)值與趨勢皆比較接近，認(rèn)為將馬斯京根法應(yīng)用于冰期是可行的。

冰蓋凍結(jié)增厚過程對洪水波傳播的影響可分為3個(gè)階段：在冰厚快速增大的第1階段，斷面綜合糙率增加，洪水波變形增大；隨后在冰厚緩慢增加的第2階段，斷面綜合糙率的降低使洪水波變形減??；最后在冰厚幾乎不發(fā)生改變但流量增大的第3階段，斷面綜合糙率沒有大的變化，流量的持續(xù)增加使洪水波變形增大。冰蓋的融化減薄過程，流量持續(xù)增加，冰厚的減小使洪水波變形增大。兩個(gè)過程的參數(shù)x均會在冰蓋凍結(jié)增厚對洪水波變形影響的第3階段有交點(diǎn)，交點(diǎn)之前的冰厚范圍，冰蓋融化減薄時(shí)的洪水波變形更大，冰下過流能力更小，交點(diǎn)之后的冰厚范圍，冰蓋融化減薄時(shí)的洪水波變形更小，冰下過流能力更大。

致謝：研究工作得到黃河水利科學(xué)研究院的大力支持，謹(jǐn)致真摯謝意！