基于馬斯京根法的多沙河流輸沙演算

楊 飛，劉博倫，江恩慧，王遠見，時志晨

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

0 引言

懸移質輸運計算是以黃河為代表的多沙河流的重要研究課題?；谒畡恿W的河流水沙數值模型是黃河干支流河道最主要的水沙計算工具，但河流水沙數值模型結構復雜、計算量大、定解條件多，不適合大流域水文模型匯總的溝道輸沙計算。Li 等［1］將人工神經網絡引入河道水沙輸運計算，與基于物理過程的水動力模型相比，減少了對地形的依賴。對于黃河這樣的大型河道洪水演進中的泥沙計算，目前還缺少快速的水文計算方法。對于一些缺乏實測地形資料的河道，當無法采用水沙數值模型計算輸沙過程時，尋找有一定精度的水文方法研究輸沙過程顯得十分必要。

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

為此，本文針對河道輸沙過程，以懸沙輸運方程和運動波方程相似性為基礎，提出一種基于馬斯京根水文學方法的泥沙輸運計算方法，并以黃河中下游干流河道水沙過程進行驗證。

1 馬斯京根法輸沙率演算

由連續方程和運動方程組成的圣維南方程組，是對河道洪水波的數學描述，連續方程積分可得河段水量平衡方程，動力方程在恒定運動波近似的前提下可得到槽蓄方程［7］，以槽蓄量與示儲流量成線性關系的假定為前提的馬斯京根流量演算目前已經在黃河干支流河道洪水演算中得到應用。

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

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

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

式中：QS，u、QS，d為河段進口和出口輸沙率，WS為河段水體內懸沙槽蓄量，Q′S為示儲輸沙率，C0～C3為演算系數，C0～C2與流量演算系數一致，SS為河道沖淤交換源項。

本文建立的流量與輸沙率的演算公式，可分段進行演算，即分段馬斯京根法。源項根據具體河道的輸沙能力與含沙量關系確定，采用不平衡輸沙形式時，源項SS為

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

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

2 北干流河道洪水演算

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

圖1 龍門—潼關段流量演算率定結果

北干流河段2018 年汛期和2019 年汛期來水量大，流量演算精度較高，確定系數分別達0.88 和0.90，演算的流量過程與實測過程基本重合。含沙量演算結果和實測結果略有偏差，確定系數分別達0.60和0.78，在不考慮河段沖淤的情況下進行含沙量演算精度是可以接受的，說明所采用的方法可以用于該河段的河道流量與含沙量的演算。

3 黃河下游水沙演算

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

圖3 小浪底至花園口河段流量驗證結果

圖4 小浪底至花園口河段含沙量驗證結果

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

圖5 花園口至夾河灘河段含沙量驗證結果

4 結論

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

本文所建立的方法高效簡單，可用該模型進行大型流域水沙模型中的溝道輸沙演算。由于以懸沙為代表的物質輸運方程具有普適性，因此所建立的方法對溶質在內的河流物質通量同樣適用，其避免了采用水動力學模型的復雜計算，在計算效率上具有優勢。