申紅彬，吳保生

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

1 研究背景

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 研究方法途徑

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 方程空時(shí)變換

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

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

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

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

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

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

對(duì)式（17）進(jìn)行變量變換，令q=VH，則可以變換為：

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

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

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

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

5 結(jié)論

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

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