楊 飛，劉博倫，江恩慧，王遠(yuǎn)見(jiàn)，時(shí)志晨

（1.黃河水利科學(xué)研究院，河南 鄭州 450003；2.水利部黃河下游河道與河口治理重點(diǎn)實(shí)驗(yàn)室，河南 鄭州 450003；3.鄭州大學(xué) 水利與交通學(xué)院，河南 鄭州 450001；4.河海大學(xué) 港口海岸與近海工程學(xué)院，江蘇 南京 210098）

0 引言

懸移質(zhì)輸運(yùn)計(jì)算是以黃河為代表的多沙河流的重要研究課題?；谒畡?dòng)力學(xué)的河流水沙數(shù)值模型是黃河干支流河道最主要的水沙計(jì)算工具，但河流水沙數(shù)值模型結(jié)構(gòu)復(fù)雜、計(jì)算量大、定解條件多，不適合大流域水文模型匯總的溝道輸沙計(jì)算。Li 等［1］將人工神經(jīng)網(wǎng)絡(luò)引入河道水沙輸運(yùn)計(jì)算，與基于物理過(guò)程的水動(dòng)力模型相比，減少了對(duì)地形的依賴(lài)。對(duì)于黃河這樣的大型河道洪水演進(jìn)中的泥沙計(jì)算，目前還缺少快速的水文計(jì)算方法。對(duì)于一些缺乏實(shí)測(cè)地形資料的河道，當(dāng)無(wú)法采用水沙數(shù)值模型計(jì)算輸沙過(guò)程時(shí)，尋找有一定精度的水文方法研究輸沙過(guò)程顯得十分必要。

應(yīng)用水文模型研究黃河流域產(chǎn)沙輸沙時(shí)，一般采用水動(dòng)力學(xué)模型或馬斯京根水文學(xué)方法進(jìn)行河網(wǎng)匯流計(jì)算［2］，多基于恒定流假定進(jìn)行平衡輸沙或不平衡輸沙計(jì)算［3－4］，與水流采用非恒定計(jì)算在理論上是矛盾的。值得注意的是，有相關(guān)學(xué)者將馬斯京根法應(yīng)用到輸沙演算過(guò)程中，Choudhury 等［5］以馬斯京根法流量演算為基礎(chǔ)，依據(jù)含沙量與流量之間的單值冪函數(shù)經(jīng)驗(yàn)關(guān)系計(jì)算含沙量，后來(lái)Sil 等［6］根據(jù)水沙單值關(guān)系，用輸沙過(guò)程推求流量過(guò)程。該方法是在平衡輸沙的框架下進(jìn)行的，嚴(yán)重依賴(lài)水沙單值關(guān)系，對(duì)于含沙量與流量之間不存在一一對(duì)應(yīng)關(guān)系的不平衡輸沙，計(jì)算則不能保證沙量守恒。

為此，本文針對(duì)河道輸沙過(guò)程，以懸沙輸運(yùn)方程和運(yùn)動(dòng)波方程相似性為基礎(chǔ)，提出一種基于馬斯京根水文學(xué)方法的泥沙輸運(yùn)計(jì)算方法，并以黃河中下游干流河道水沙過(guò)程進(jìn)行驗(yàn)證。

1 馬斯京根法輸沙率演算

由連續(xù)方程和運(yùn)動(dòng)方程組成的圣維南方程組，是對(duì)河道洪水波的數(shù)學(xué)描述，連續(xù)方程積分可得河段水量平衡方程，動(dòng)力方程在恒定運(yùn)動(dòng)波近似的前提下可得到槽蓄方程［7］，以槽蓄量與示儲(chǔ)流量成線性關(guān)系的假定為前提的馬斯京根流量演算目前已經(jīng)在黃河干支流河道洪水演算中得到應(yīng)用。

引入示儲(chǔ)流量與河段槽蓄量成單一線性關(guān)系的假設(shè)，是馬斯京根流量演算法的基本出發(fā)點(diǎn)。馬斯京根流量演算法由河段水量平衡方程［式（1）］和槽蓄方程［式（2）］聯(lián)立求解，獲得河段流量演算公式［式（3）］：

式中：Qu、Qd為河段進(jìn)口和出口流量，上標(biāo)n和n＋1 分別表示t時(shí)刻和t＋Δt時(shí)刻，Δt為計(jì)算時(shí)段長(zhǎng)，W為河段槽蓄量，Q′為示儲(chǔ)流量，C0、C1、C2均為流量演算系數(shù)，K為蓄量常數(shù)（相當(dāng)于洪水波在河段中的傳播時(shí)間），x為流量比重因數(shù)（主要與洪水波的坦化變形程度有關(guān)）。

不考慮慣性項(xiàng)時(shí)，運(yùn)動(dòng)波方程的空間偏心四點(diǎn)離散求解方法具有馬斯京根法的形式，當(dāng)運(yùn)動(dòng)波方程有限差分?jǐn)?shù)值解的誤差剛好等于擴(kuò)散波的物理擴(kuò)散項(xiàng)時(shí)，該數(shù)值解即為擴(kuò)散波方程的精確解［8］，這為該方法的擴(kuò)散波洪水演算提供了一定的理論支撐。運(yùn)動(dòng)波方程的波速為不考慮沖淤與擴(kuò)散時(shí)一維懸沙輸運(yùn)方程為公式形式與運(yùn)動(dòng)波方程完全相同，對(duì)應(yīng)的波速為與洪水運(yùn)動(dòng)波的波速略有差異，當(dāng)忽略?xún)烧卟町悤r(shí)，兩者的演進(jìn)特征一致。因此，理論上適用于洪水運(yùn)動(dòng)波和擴(kuò)散波演算的馬斯京根法同樣適用于泥沙輸運(yùn)演算。假定不論在輸沙率增大階段還是減小階段，河段內(nèi)輸沙率是沿程線性變化的，當(dāng)示儲(chǔ)流量與河段槽蓄量成單一關(guān)系時(shí)，采用與Q′同比重的輸沙率Q′S，與該河段水體懸沙槽蓄量WS亦能成單一線性關(guān)系，稱(chēng)式（5）為馬斯京根法的輸沙槽蓄關(guān)系式或輸沙槽蓄方程，這里輸沙的蓄量常數(shù)假定和比重因數(shù)與水流一致。河段內(nèi)的泥沙同樣滿(mǎn)足質(zhì)量平衡，單位時(shí)間內(nèi)的泥沙質(zhì)量變化量等于進(jìn)口與出口輸沙率差值加上河道沖淤源項(xiàng)，因此有沙量平衡方程［式（6）］。式（6）和式（5）聯(lián)立求解獲得馬斯京根法的輸沙率演算公式［式（7）］，本文的輸沙率演算思想與Singh 等［9］所建立潰壩洪水演進(jìn)模型輸沙計(jì)算方法一致。

式中：QS，u、QS，d為河段進(jìn)口和出口輸沙率，WS為河段水體內(nèi)懸沙槽蓄量，Q′S為示儲(chǔ)輸沙率，C0～C3為演算系數(shù)，C0～C2與流量演算系數(shù)一致，SS為河道沖淤交換源項(xiàng)。

本文建立的流量與輸沙率的演算公式，可分段進(jìn)行演算，即分段馬斯京根法。源項(xiàng)根據(jù)具體河道的輸沙能力與含沙量關(guān)系確定，采用不平衡輸沙形式時(shí)，源項(xiàng)SS為

式中：S、S?分別為河段平均懸移質(zhì)泥沙濃度和挾沙力，ω為泥沙顆粒的沉速，α為懸移質(zhì)泥沙的恢復(fù)飽和系數(shù)，h為河道斷面平均水深，L為河道長(zhǎng)度，A為河道斷面過(guò)流面積。

也可根據(jù)實(shí)際河段情況采用經(jīng)驗(yàn)公式計(jì)算源項(xiàng)。當(dāng)不考慮源項(xiàng)時(shí)，輸沙率與流量的演算公式形式一致。本方法中的水沙傳播時(shí)間對(duì)應(yīng)了水流與泥沙的運(yùn)動(dòng)速度，本文在計(jì)算泥沙輸運(yùn)時(shí)，采用了與水流相同的傳播時(shí)間，這意味著本方法假定兩者同步傳播。當(dāng)兩者傳播速度不一致時(shí)，導(dǎo)致沙峰滯后于洪峰，需要選擇泥沙傳播時(shí)間進(jìn)行計(jì)算。

2 北干流河道洪水演算

對(duì)2018 年汛期和2019 年汛期北干流河段的水沙過(guò)程進(jìn)行演算，選取河曲、府谷、吳堡、龍門(mén)、潼關(guān)等5個(gè)水文站為節(jié)點(diǎn)將河曲到潼關(guān)河段分為4 段分別演算，以每段上游節(jié)點(diǎn)實(shí)測(cè)水沙過(guò)程為輸入條件，計(jì)算下游節(jié)點(diǎn)的水沙過(guò)程，其中龍門(mén)—潼關(guān)段考慮支流匯流過(guò)程。北干流河段內(nèi)有天橋水電站，庫(kù)容較小，對(duì)水流有一定的調(diào)節(jié)能力，能夠改變中小流量的水沙過(guò)程，本文沒(méi)有考慮這部分影響。實(shí)測(cè)流量過(guò)程用于模型參數(shù)的率定，實(shí)測(cè)輸沙過(guò)程用于模擬驗(yàn)證。龍門(mén)—潼關(guān)段的流量和含沙量過(guò)程計(jì)算結(jié)果如圖1 和圖2 所示。

圖1 龍門(mén)—潼關(guān)段流量演算率定結(jié)果

北干流河段2018 年汛期和2019 年汛期來(lái)水量大，流量演算精度較高，確定系數(shù)分別達(dá)0.88 和0.90，演算的流量過(guò)程與實(shí)測(cè)過(guò)程基本重合。含沙量演算結(jié)果和實(shí)測(cè)結(jié)果略有偏差，確定系數(shù)分別達(dá)0.60和0.78，在不考慮河段沖淤的情況下進(jìn)行含沙量演算精度是可以接受的，說(shuō)明所采用的方法可以用于該河段的河道流量與含沙量的演算。

3 黃河下游水沙演算

將所建立的模型用于黃河小浪底至利津河段的水沙演算。在小浪底至花園口河段，考慮伊洛河、沁河入?yún)R，以黑石關(guān)、武陟作為支流入口，根據(jù)洪水傳播時(shí)間設(shè)定計(jì)算時(shí)段長(zhǎng)為13 h，采用2006 年6 月1 日至2017年汛期實(shí)測(cè)流量用于率定參數(shù)，最小二乘法率定結(jié)果的確定系數(shù)為0.98，對(duì)2018 年汛期至2021 年汛期流量過(guò)程驗(yàn)證，整個(gè)時(shí)段驗(yàn)證結(jié)果的確定系數(shù)為0.98，圖3 為2018 年汛期和2019 年汛期流量過(guò)程的驗(yàn)證結(jié)果，可以看出，馬斯京根法對(duì)該河段流量演算的精度較高，演算流量過(guò)程與實(shí)測(cè)過(guò)程基本一致。對(duì)2006 年汛期至2021 年汛期的輸沙率進(jìn)行驗(yàn)證，確定系數(shù)為0.57，圖4 為2018 年汛期和2019 年汛期含沙量過(guò)程的驗(yàn)證結(jié)果，演算的含沙量過(guò)程整體上與實(shí)際過(guò)程擬合較好，總體上看模型可以用于小浪底至花園口河段的河道流量與含沙量演算。小浪底至花園口河段為典型的游蕩型河段，洪水期河道沖淤幅度大，輸沙調(diào)整明顯且與含沙量關(guān)系密切［10］，導(dǎo)致含沙量較大時(shí)演算過(guò)程與實(shí)測(cè)過(guò)程偏差較大，為了提高輸沙率的演算精度，需進(jìn)一步在輸沙率計(jì)算中考慮泥沙沖淤源項(xiàng)。

圖3 小浪底至花園口河段流量驗(yàn)證結(jié)果

圖4 小浪底至花園口河段含沙量驗(yàn)證結(jié)果

對(duì)于花園口至利津河段，選取花園口、夾河灘、高村、孫口、艾山、濼口、利津等7 個(gè)水文站2018 年汛期和2019 年汛期的水沙數(shù)據(jù)進(jìn)行分析驗(yàn)證。以7 個(gè)水文站為節(jié)點(diǎn)將花園口至利津河段分為6 段分別進(jìn)行洪水演算，不考慮區(qū)間內(nèi)沖淤、引水、支流入?yún)R對(duì)水量和沙量的影響，演算下游節(jié)點(diǎn)的水沙過(guò)程。6 個(gè)節(jié)點(diǎn)的流量演算精度較高，確定系數(shù)均在0.95 以上，這里不再展示。2018 年和2019 年夾河灘站汛期含沙量驗(yàn)證結(jié)果如圖5 所示，確定系數(shù)分別為0.91 和0.84，其他5個(gè)節(jié)點(diǎn)含沙量驗(yàn)證結(jié)果的確定系數(shù)均在0.90 以上?？梢钥闯龈鞫蔚暮沉垦菟阒蹬c實(shí)測(cè)值非常接近，采用馬斯京根方法有很高的精度。受泥沙沿程沖淤變化等影響，泥沙輸運(yùn)的模擬精度略低于流量演算的。

圖5 花園口至夾河灘河段含沙量驗(yàn)證結(jié)果

4 結(jié)論

一維懸沙輸運(yùn)方程和運(yùn)動(dòng)波方程的形式一致且波速接近，以運(yùn)動(dòng)波方程為理論基礎(chǔ)的馬斯京根流量演算方法同樣適用于懸沙輸運(yùn)演算。據(jù)此，本文參照馬斯京根流量演算方法，建立了馬斯京根法輸沙率演算方法。將該方法用于黃河中下游干流水沙過(guò)程的演算，結(jié)果表明，基于馬斯京根法的輸沙演算方法能保證一定精度并且計(jì)算過(guò)程十分簡(jiǎn)單，可應(yīng)用在以懸移輸沙為主的多沙河流水沙演算。水文預(yù)報(bào)中當(dāng)洪水流量過(guò)程采用馬斯京根法進(jìn)行演算時(shí)，該方法參數(shù)能夠直接用于洪水輸沙預(yù)報(bào)。

本文所建立的方法高效簡(jiǎn)單，可用該模型進(jìn)行大型流域水沙模型中的溝道輸沙演算。由于以懸沙為代表的物質(zhì)輸運(yùn)方程具有普適性，因此所建立的方法對(duì)溶質(zhì)在內(nèi)的河流物質(zhì)通量同樣適用，其避免了采用水動(dòng)力學(xué)模型的復(fù)雜計(jì)算，在計(jì)算效率上具有優(yōu)勢(shì)。