河流泥沙水文學(xué)模型邊界條件參數(shù)化方法探討

申紅彬，吳保生

（1.華北水利水電大學(xué)，河南 鄭州 450045；2.北京師范大學(xué) 水科學(xué)研究院 城市水循環(huán)與海綿城市技術(shù)北京市重點實驗室，北京 100875；3.清華大學(xué) 水沙科學(xué)與水利水電工程國家重點實驗室，北京 100084）

1 研究背景

河流泥沙數(shù)學(xué)模型主要包括：水文學(xué)模型、水動力學(xué)模型和水文水動力學(xué)模型［1］。其中，水文學(xué)模型主要以分析河道實測水文泥沙資料為基礎(chǔ)，建立經(jīng)驗關(guān)系或概念模型。水文學(xué)模型往往因形式簡單、方便實用，因而在生產(chǎn)實踐中得到了較為廣泛的應(yīng)用?，F(xiàn)有水文學(xué)模型受參數(shù)不確定性影響，經(jīng)常會面臨以下應(yīng)用問題：（1）延用性問題，即原有河道模型參數(shù)是否適用于變化環(huán)境；（2）移植性問題，即某河道模型參數(shù)能否推廣應(yīng)用于其它河道。水文學(xué)模型方法雖然不像傳統(tǒng)經(jīng)典數(shù)學(xué)物理方程那樣具有嚴密的物理概念與數(shù)學(xué)推導(dǎo)，但作為對客觀世界的另一種描述方式，也應(yīng)具有一定的物理基礎(chǔ)［2-3］。因此，要想發(fā)展完善河流泥沙水文學(xué)模型就需要解決參數(shù)的確定方法以及方程的物理意義問題。

現(xiàn)有河流泥沙水文學(xué)模型主要以考慮上站來水含沙量的輸沙冪律函數(shù)公式：（S為河段出口含沙量，Q、S0分別為河段流量、進口來水含沙量，K為輸沙系數(shù)，a、b為輸沙指數(shù)）作為基礎(chǔ)方程，系數(shù)K與指數(shù)a、b主要用以反映河流邊界條件的影響。其中，系數(shù)K主要與河床前期累積沖淤量有關(guān)［4］，指數(shù)a、b主要與河槽斷面形態(tài)、沿程距離等幾何因素有關(guān)［5-7］。吳保生等［8］在對黃河內(nèi)蒙古河段淤積計算過程中發(fā)現(xiàn)，考慮河床沖淤演變過程中系數(shù)K值變化可有效提高河槽淤積量的計算精度，并認為指數(shù)a、b值變化實質(zhì)反映了輸沙能力在空間上隨河床邊界條件的變化，系數(shù)K值變化實質(zhì)反映了輸沙能力在時間上隨河床邊界條件調(diào)整的變化。王彥君等［9］進一步考慮了河道年內(nèi)汛期、非汛期輸沙指數(shù)a、b的差異。姚文藝等［5］通過模型試驗直觀反映了黃河下游河道不同萎縮模式（“集中淤槽”與“灘槽并淤”）下輸沙指數(shù)a、b的調(diào)整變化情況。鑒于系數(shù)K與指數(shù)a、b對輸沙與沖淤計算精度的重要影響，不少學(xué)者［10-12］試圖通過統(tǒng)計回歸得到不同河段模型參數(shù)（系數(shù)K、指數(shù)a、b）值。不過，這種分析方法偏于經(jīng)驗，得到不同河段、不同時期內(nèi)的系數(shù)K、指數(shù)a、b值變化范圍較大，存在一定的不確定性。

本文在分析河流泥沙水文學(xué)模型經(jīng)驗關(guān)系方程基本內(nèi)涵的基礎(chǔ)上，結(jié)合對現(xiàn)有模型參數(shù)研究成果的分析比較，提出根據(jù)河道邊界條件確定模型參數(shù)的方法，概括總結(jié)河道邊界條件參數(shù)化的研究途徑，并從河道邊界條件參數(shù)化角度出發(fā)，通過方程空間—時間變量變換，分析探討河流泥沙水文學(xué)模型方程的物理意義。

2 參數(shù)內(nèi)涵分析

作為河流泥沙水文學(xué)模型的基礎(chǔ)方程，考慮上站來水含沙量的輸沙冪律函數(shù)公式可以表示為［4，12-13］：

式中：S為河段出口斷面含沙量，kg/m3；Q為流量，m3/s；S0為河段進口斷面來水含沙量，kg/m3；K為輸沙系數(shù)，與河床前期累計沖淤量有關(guān)；a、b為輸沙指數(shù)，統(tǒng)計表明指數(shù)a、b滿足：a+b≈1.0［14-16］。

分析河流泥沙水文學(xué)模型基礎(chǔ)方程式（1），本身蘊含了以下基本內(nèi)涵：

（1）體現(xiàn)了不平衡輸沙的理念。當(dāng)來水含沙量大于出口含沙量時，河道處于淤積狀態(tài)；當(dāng)來水含沙量小于出口含沙量時，河道處于沖刷狀態(tài)；當(dāng)來水含沙量等于出口含沙量時，河道處于輸沙平衡狀態(tài)。當(dāng)河道處于輸沙平衡狀態(tài)時，來水含沙量與出口含沙量應(yīng)相等于水流挾沙力，則有：

式中S*為河道水流挾沙力，也稱臨界含沙量。

因指數(shù)a、b滿足統(tǒng)計關(guān)系：a+b≈1.0，相應(yīng)式（2）中流量指數(shù)也可表示為：b/（1-a）≈1.0，這表明式（2）流量指數(shù)b/（1-a）具有統(tǒng)計平均意義上的穩(wěn)定性。

結(jié)合式（2），暫時忽略河段沿程流量變化，則基礎(chǔ)方程式（1）可表示為其它形式：

式中：S/S*為河道沿程斷面含沙量與臨界含沙量的比值；S0/S*為河道進口斷面含沙量與臨界含沙量的比值；SDR為河段排沙比；S0/Qb/（1-a）為來沙系數(shù)。

式（3）左側(cè)公式表明，考慮上站來水含沙量的水沙關(guān)系模型實質(zhì)反映了河道斷面含沙量與臨界含沙量比值的沿程變化情況。根據(jù)不平衡輸沙理論，可以判斷指數(shù)a值應(yīng)在［0，1］之間，當(dāng)河道進出口斷面距離趨于0時指數(shù)a值應(yīng)等于1.0，當(dāng)進出口斷面距離無窮大時指數(shù)a值趨于0［17］。

（2）反映了相對空間的概念。根據(jù)模型方程，對于同一出口斷面，當(dāng)上游來水含沙量取不同距離進口斷面時，模型參數(shù)會有所區(qū)別；同理，對于同一進口斷面，下游不同距離出口斷面的輸沙參數(shù)也有所不同。假設(shè)對于同一河道可以劃分為0-1、1-2、…、（n-1）-n共計n個河段，各河段模型參數(shù)相同均為K、a、b，則應(yīng)用模型方程依次遞推，可以得到以下關(guān)系：

式中：n為河道劃分河段數(shù)；i為迭代計算序號。

從式（4）可以看出，對于模型方程式（1），同一進口斷面下游不同距離出口斷面的模型系數(shù)、指數(shù)是不同的。不過，對于下游每一個出口斷面，模型方程系數(shù)、指數(shù)均具有如下關(guān)系：

式（5）表明，對于下游每一個出口斷面，河道水流挾沙能力方程式（2）的系數(shù)、指數(shù)是一致的。

（3）賦予了反映河流自動調(diào)整作用的功能。通過建立系數(shù)K與河槽累計沖淤量的變化關(guān)系，當(dāng)河道發(fā)生累計淤積時，K隨累計淤積量的增加而增大；反之當(dāng)河道發(fā)生累計沖刷時，K隨累計沖刷量的增加而減小，這實質(zhì)體現(xiàn)了河槽邊界條件調(diào)整對輸沙能力的影響［8，10，18］：

式中：K0為輸沙系數(shù)基礎(chǔ)值；λ為系數(shù)；為前期河床沖淤參數(shù)，億t；ΔWs為計算時段沖淤量，億t；D為平均沖淤參數(shù)，億t。

綜合以上內(nèi)涵分析可以看出，對于河流泥沙水文學(xué)模型基礎(chǔ)方程：，模型系數(shù)、指數(shù)更多受到河道邊界條件的影響。所謂河流泥沙水文學(xué)模型邊界條件參數(shù)化方法，就是以河流泥沙水文學(xué)模型基礎(chǔ)方程為基礎(chǔ)，從河道邊界條件（如河段長度、寬度、斷面形態(tài)、河床比降等）參數(shù)化角度出發(fā)，根據(jù)河道邊界條件確定模型的系數(shù)、指數(shù)參數(shù)，可表達為如下函數(shù)形式：

式中L、B、H、J分別為反映河道邊界條件的長度、寬度、水深及比降。

3 研究方法途徑

3.1 統(tǒng)計分析方法對式（1）兩側(cè)取對數(shù)，可以得到：

基于式（8），根據(jù)河道沿程不同斷面實測水沙資料，通過多元線性回歸，可以得到不同河段模型參數(shù)（系數(shù)K與指數(shù)a、b）值，進而分析模型參數(shù)的變化規(guī)律，這就是統(tǒng)計分析方法。

圖1為不少研究者基于黃河沿程水文站（圖1（a））實測水沙資料統(tǒng)計回歸得到的模型參數(shù)（系數(shù)K與指數(shù)a、b）值變化情況（圖1（b）），各水文站河段進口斷面為上游鄰近水文站，具體參見文獻［12，19-20］）。不過，因不同河段邊界條件相差甚遠，且系數(shù)K值又與河床前期累積沖淤條件有關(guān)，河段選取缺乏統(tǒng)一的標(biāo)準(zhǔn)，導(dǎo)致不同河段模型參數(shù)（系數(shù)K與指數(shù)a、b）變化情況復(fù)雜。其中，指數(shù)和（a+b）的統(tǒng)計平均值約為1.1，接近于1.0。

根據(jù)模型方程內(nèi)涵（2）的分析，河段選取宜先對同一河段固定同一進口斷面，分別選取下游不同距離出口斷面，基于實測水沙資料統(tǒng)計回歸模型參數(shù)。圖1（c）為根據(jù)式（3）右側(cè)排沙比公式，以黃河下游花園口站為同一進口斷面，分別選取下游不同距離出口斷面：高村、艾山、利津，簡單統(tǒng)計回歸河段排沙比與花園口來沙系數(shù)之間關(guān)系的變化情況（圖中：SDRgc、SDRas、SDRlj分別為花園口-高村、花園口-艾山、花園口-利津河段排沙比；Sh/Q為花園口站來沙系數(shù)，b/（1-a）≈1.0）。從中可以看出，隨著河段距離的增大，模型系數(shù)K、指數(shù)a的值逐步減小。這為利用統(tǒng)計分析方法直接分析模型參數(shù)隨距離的變化規(guī)律提供了方便。后再對不同河段比較分析距離因子對輸沙參數(shù)的影響規(guī)律。費祥俊等［17］基于黃河下游河道實測水沙資料，通過統(tǒng)計回歸，得到艾山至利津、三黑小（三門峽、黑石關(guān)、小董）至艾山、三黑小至利津河段排沙比公式（3）指數(shù)（a-1）值依次為-0.03、-0.50、-0.53，相應(yīng)指數(shù)a值依次為：0.97、0.50、0.47，顯示出河段距離越大指數(shù)a值越小的趨勢。

鑒于天然河道邊界特征沿程變化多樣，可基于渠道觀測或水槽試驗資料，通過統(tǒng)計分析，研究考慮上站來水含沙量的河道水沙關(guān)系模型參數(shù)（系數(shù)K與指數(shù)a、b）值的變化規(guī)律。例如，圖2為長江科學(xué)研究院水槽沖刷試驗測得的含沙量沿程恢復(fù)情況（圖2（a））［21］，固定選取某一斷面作為進口斷面，基于式（3）左側(cè)公式可以率定出下游不同距離出口斷面指數(shù)a值（圖2（b））?？梢钥闯?，隨著河段出口斷面距離的增大，輸沙指數(shù)a值呈現(xiàn)出明顯的衰減趨勢。不過，目前對于考慮上站來水含沙量的輸沙冪律函數(shù)公式，多用于天然河道特別是多沙河流泥沙輸移的模擬，較少用于描述渠道或水槽內(nèi)泥沙的輸移。

圖1 黃河不同河段模型參數(shù)統(tǒng)計回歸值變化情況

圖2 基于水槽試驗資料的模型參數(shù)統(tǒng)計分析方法

3.2 理論比較途徑河流泥沙三類模型（水文學(xué)模型、水動力學(xué)模型、水文水動力學(xué)模型）作為對河流泥沙輸移同一物理現(xiàn)象的不同描述，在表達形式上雖存在區(qū)別但也存在聯(lián)系。通過對3類模型開展理論比較分析，可用于探求河流泥沙水文模型參數(shù)的計算方法。

在河流泥沙水文水動力學(xué)模型中，水流挾沙力方程一般采用基于水流功率理論的簡化形式：QS*=QS*=A（QJ）m（A為系數(shù)，m為指數(shù)，J為比降）［22］。其中，指數(shù)m值一般認為在2.0左右，對于不同斷面形態(tài)河槽指數(shù)m值存在一定的區(qū)別。比較水流挾沙力方程式（2），兩者方程結(jié)構(gòu)形式存在相似性。因此，對于水流挾沙力方程式（2）中的流量指數(shù)b/（1-a）值，應(yīng)等于（m-1）約為1.0左右，這與實測資料統(tǒng)計平均結(jié)果相符，而對于系數(shù)值K1/（1-a），判斷應(yīng)與河道比降有關(guān)。

在水動力學(xué)模型中，一維恒定流不平衡輸沙微分方程為：

式中：x為沿程距離；ω為泥沙沉速；q為單寬流量；，sb*為河底水流挾沙力，也稱為泥沙恢復(fù)飽和系數(shù)。

根據(jù)河流泥沙水文學(xué)模型基礎(chǔ)方程內(nèi)涵（1）、（2）的討論，并參考結(jié)合圖2，指數(shù)a應(yīng)是一個與輸沙距離有關(guān)的變量，可以近似采用指數(shù)函數(shù)進行描述：a=e-κx，相應(yīng)輸沙微分方程可表示為［23］：

式中κ為變化速率參數(shù)。

筆者曾對式（9）、式（10）進行比較［23］，通過對Sln（S/S*）在S*附近進行一階Taylor級數(shù)展開，表明Sln（S/S*）與（S-S*）近似相等：

因此，式（10）變化速率參數(shù)κ=α*ω/q，相應(yīng)指數(shù)a計算表達式為：

根據(jù)式（12），相同流量、比降情況下，不同斷面形態(tài)河槽下游不同距離斷面泥沙水文模型指數(shù)a及含沙量沿程變化示意情況如圖3所示。

圖3 相同流量、比降下不同河寬河槽泥沙水文模型指數(shù)a與含沙量沿程變化示意圖

基于式（12），輸沙微分方程式（10）可進一步表示為：

至于輸沙微分方程式（9）與式（10）、式（13）何者描述形式更為合理，尚需后期進一步的研究。

4 方程空時變換

對于一維恒定流，對式（9）、式（13）進行變量變換，令dx=Vdt、q=VH（其中，V為均勻流平均流速，H為平均水深），則二式可以分別變換為：

式（14）、式（15）可以統(tǒng)一表示為如下形式：

式中：y為系統(tǒng)特征變量；t為時間；ye為系統(tǒng)特征變量平衡值；β為變化速率。

因此，從“系統(tǒng)”角度來看，式（14）、式（15）均為系統(tǒng)滯后響應(yīng)的微分方程，系統(tǒng)特征變量y調(diào)整變化速率與該特征變量當(dāng)前狀態(tài)y和平衡狀態(tài)ye之間的差值成正比［24-26］。從斷面質(zhì)點系統(tǒng)跟蹤角度來看，式（14）、式（15）描述的是河道水沙質(zhì)點系統(tǒng)在進口水沙初始條件下，經(jīng)過一定時間（對應(yīng)一定運移距離）后，自身水沙組合情況調(diào)整并趨向于平衡的過程（見圖4），不同之處在于兩者對于描述系統(tǒng)水沙組合的特征變量y的選擇不同，分別為（S/S*）與ln（S/S*）。

圖4 河道斷面水沙質(zhì)點系統(tǒng)滯后響應(yīng)過程

在一維恒定均勻流情況下，微分方程式（9）的解析解為：

對式（17）進行變量變換，令q=VH，則可以變換為：

式中：t為時間，t=x/V；T為泥沙沉降時間，T=H/ω。

從式（18）可以看出，以往微分方程式（9）對應(yīng)的含沙量關(guān)于空間變化的解析解，經(jīng)過變量變換，也可以表示為關(guān)于時間變化的滯后響應(yīng)解的形式。

同理，基于式（12），對于河流泥沙水文學(xué)模型方程式（3），經(jīng)過變量變換也可以表示為：

從式（19）可以看出，原用于確定指數(shù)a的河道地貌變量：河長x、河寬B，經(jīng)變量變換為：河長x、水深H后，再分別與流速V、泥沙沉速ω相除，最終轉(zhuǎn)化為系統(tǒng)自身調(diào)整時間t與斷面泥沙沉降時間T，這可以從另一角度理解河流泥沙水文模型基礎(chǔ)方程的物理意義。

5 結(jié)論

（2）河流泥沙水文學(xué)模型邊界條件參數(shù)化研究方法途徑主要包括：統(tǒng)計分析方法與理論比較途徑。統(tǒng)計分析方法宜先對同一河段選用同一進口斷面，分別選取下游不同距離出口斷面，基于實測水沙資料統(tǒng)計回歸模型參數(shù)，進而研究參數(shù)隨距離變化規(guī)律，后再對不同河段比較分析距離因子對輸沙參數(shù)的影響規(guī)律。鑒于天然河道邊界特征沿程變化多樣，可先選用渠道觀測或水槽試驗資料，統(tǒng)計分析模型參數(shù)的變化規(guī)律。理論比較途徑主要基于對水文學(xué)模型、水動力學(xué)模型以及水文水動力學(xué)模型三種模型方程的理論比較分析，探求河流泥沙水文學(xué)模型參數(shù)的計算方法。

（3）從斷面質(zhì)點系統(tǒng)跟蹤角度來看，描述河流泥沙運動的水文學(xué)與水動力學(xué)模型關(guān)于空間變化的數(shù)學(xué)方程形式，經(jīng)過變量變換均可表示為關(guān)于時間變化的滯后響應(yīng)模型方程形式，兩者區(qū)別在于對描述系統(tǒng)水沙組合的特征變量選擇不同。其中，對于河流泥沙水文模型邊界條件參數(shù)化的最終結(jié)果轉(zhuǎn)化為了水沙質(zhì)點系統(tǒng)的自身調(diào)整時間與斷面泥沙沉降時間兩個變量，從另一角度提供了對河流泥沙水文學(xué)模型基礎(chǔ)方程物理意義的理解。