周美蓉，夏軍強，鄧珊珊，劉 鑫

（武漢大學(xué) 水資源與水電工程科學(xué)國家重點實驗室，湖北 武漢 430072）

1 研究背景

水流挾沙力通常指在一定的水流及泥沙綜合條件下，水流能夠攜帶懸移質(zhì)中床沙質(zhì)的臨界含沙量［1］。水流挾沙力公式及其參數(shù)選取的合理性，直接影響到河床沖淤變形的計算精度。目前國內(nèi)外已有眾多學(xué)者對水流挾沙力公式進(jìn)行了深入研究，提出各類經(jīng)驗或半經(jīng)驗半理論的計算公式［2-5］。在少沙河流上，如長江中下游，通常采用張瑞瑾公式計算水流挾沙力，其結(jié)構(gòu)形式簡明，實際應(yīng)用方便［5-8］。該公式一般可寫成如下形式：

式中：k為包含量綱的系數(shù)，kg/m3；m為指數(shù)；ωm為非均勻懸沙的平均沉速，m/s。此處定義U3/ (ghωm)為水沙綜合參數(shù)，并用符號C表示。

然而應(yīng)用張瑞瑾挾沙力公式的實際困難在于公式中系數(shù)k和指數(shù)m的取值不易確定，兩者在不同水沙條件下存在不同取值，需要根據(jù)實測資料確定。目前確定各類挾沙力公式中參數(shù)的方法主要包括兩類：通過不斷調(diào)試系數(shù)和指數(shù)，使得計算的含沙量過程與實測值能較好地符合［7-9］；通過大量實測挾沙力資料，直接率定出公式中的系數(shù)與指數(shù)［10-11］。目前采用第二類研究方法的成果較少，如武漢水利電力學(xué)院水流挾沙力研究組［12］曾采用1956—1958年長江中下游資料（僅有1 組含沙量小于0.1 kg/m3）對張瑞瑾挾沙力公式進(jìn)行率定（k≈0.0530，m≈1.54），而三峽工程運用后，懸沙含量急劇減小，且以小于0.1 kg/m3為主，故有必要增加超低含沙量數(shù)據(jù)，使率定的參數(shù)更適合三峽工程運行后的情況；方波［10］則收集了長江下游大通-鎮(zhèn)江河道大量水文泥沙實測資料，采用|Δh/Δt|≤0.1 m/d 作為判斷沖淤相對平衡的指標(biāo)，選取實測水流挾沙力數(shù)據(jù)，由此率定出該河段張瑞瑾挾沙力公式中的參數(shù)（k=0.068 kg/m3，m=1.457）；Xia 等［11］根據(jù)黃河下游8 個水文站的實測挾沙力資料，驗證了張紅武挾沙力公式中參數(shù)取值的準(zhǔn)確性（k=2.5 kg/m3，m=0.62），該公式多用于高含沙河流。

但不同研究者通過上述兩種方法確定的參數(shù)取值差異較大。余明輝等［13］采用張瑞瑾挾沙力公式時，k和m分別取為0.25 kg/m3和0.75；在郭慶超［14］和Yuan 等［8］的研究中，系數(shù)k為0.07 kg/m3，指數(shù)m為1.14；在Zhou 等［7］的模型中，k和m則分別為0.15 kg/m3和1.0。主要是由于不同研究河段與時段內(nèi)的水流、泥沙等條件不同，導(dǎo)致參數(shù)的取值亦不同。因此，對某一長河段或同一河段的長時段模擬，若將挾沙力公式的參數(shù)取為常數(shù)，在水力泥沙要素變化較大時，數(shù)學(xué)模型難以準(zhǔn)確地模擬水沙輸移過程。故有必要根據(jù)實際的水沙條件，調(diào)整挾沙力公式的參數(shù)取值，從而提高數(shù)學(xué)模型應(yīng)用于天然河流的適應(yīng)能力。實際上，張瑞瑾等［1］的研究成果已表明，挾沙力公式的系數(shù)和指數(shù)隨水沙綜合參數(shù)C而變化，并給出了參數(shù)k，m與C之間的關(guān)系（圖1）。圖中曲線由許多準(zhǔn)平衡河流實測資料和水槽試驗數(shù)據(jù)率定得到（原武漢水利電力學(xué)院玻璃水槽、南京水利實驗處鋼板水槽、長江、黃河等）。由圖1可知，當(dāng)水沙綜合參數(shù)C小于4.0 時，指數(shù)m為常數(shù)，約等于1.5，而當(dāng)C大于4.0 時，m隨其增加而減小；系數(shù)k則一直隨著C的增加而增加，但在C小于10.0 且含沙量小于0.1 kg/m3的范圍內(nèi)，k的變化過程未給出。然而三峽工程運用后水沙綜合參數(shù)C值基本小于10.0，且多年平均含沙量由工程運用前的1.0 kg/m3減小到0.1 kg/m3左右。故需重新率定低含沙量及小水沙綜合參數(shù)條件下張瑞瑾公式中的系數(shù)k和指數(shù)m，以適用于近期長江中游懸沙輸移計算。

圖1 挾沙力公式中參數(shù)k、m 與水沙綜合參數(shù)C 之間的關(guān)系［1］

基于上述分析，本研究收集了大量實測資料，開展了長江中游河段水流挾沙力的研究，確定了低含沙量（即含沙量小于10 kg/m3）及小水沙綜合參數(shù)情況下張瑞瑾挾沙力公式中參數(shù)的取值方法，并應(yīng)用于長江中游荊江河段的水沙輸移計算。

2 低含沙量水流挾沙力資料整理

在本次分析中，收集了長江中游河段的三組實測資料，用于研究水流挾沙力與相應(yīng)水沙條件之間的關(guān)系。各測站水力泥沙要素變化范圍，如表1所示。

第一組：三峽工程運用前，長江中游宜昌、陳家灣、沙市、新廠、監(jiān)利、洪水港、螺山、漢口、青山（南）及青山（北）10 個水文或水位站1956—1958年各測次的水力泥沙要素，排除數(shù)據(jù)缺失或偏差較大的23 測次，共計80 組［15］。這些實測數(shù)據(jù)包括懸移質(zhì)中床沙質(zhì)含沙量S′（扣除沖瀉質(zhì)部分，0.0026 ～0.4110 kg/m3），流量Q（3720 ～70 500 m3/s），流速U（0.75 ～2.79 m/s），水深h（3.1 ～17.9 m），平均沉降速度ω′m（扣除沖瀉質(zhì)）等。該組數(shù)據(jù)測量時盡量保證了河床為基本沖淤平衡狀態(tài)，并經(jīng)過嚴(yán)格的審查和挑選，實測含沙量大小總體上可近似等于水流挾沙力。

第二組：三峽工程運用前，長江中游新廠站1982—1985年數(shù)據(jù)共40 測次，包括懸移質(zhì)含沙量S（未扣除沖瀉質(zhì)部分，0.093 ～3.760 kg/m3），流速U（0.77 ～2.10 m/s），水深h（3.6 ～13.3 m），懸沙級配ΔPsk，床沙級配ΔPbk等［16］。這組數(shù)據(jù)也在河床基本沖淤平衡時測得。

第三組：三峽工程運用后，長江中游荊江段枝城、沙市和監(jiān)利水文站的實測資料共453 測次，包括懸移質(zhì)含沙量S（未扣除沖瀉質(zhì)部分，0.002 ～0.812 kg/m3），流速U（0.43 ～2.43 m/s），水深h（4.5 ～17.1 m），水位Z（23.91 ～46.18 m），懸沙級配ΔPsk，床沙級配ΔPbk等。三峽工程的運用使壩下游河段經(jīng)歷持續(xù)的河床沖刷，但存在某些特定時段，河床處于相對平衡狀態(tài)，此時的含沙量可近似等于挾沙力。本研究判定河床是否處于相對沖淤平衡狀態(tài)的標(biāo)準(zhǔn)為平均床面高程的變化率小于0.01 m/d。通過水位減去水深得到平均床面高程，則為平均河床高程的變化量，即為平均床面高程變化率。

需說明，較短時段內(nèi)的河床沖淤狀態(tài)不能反映其真實情況，故武漢水利電力學(xué)院水流挾沙力研究組［15］指出：需要進(jìn)一步計算各測站相鄰三天的沖淤變化，來確保該斷面處于相對平衡狀態(tài)。而在本研究中，水沙要素的測量時間不連續(xù)，無法采取上述方法進(jìn)行沖淤平衡的數(shù)據(jù)審查。但此處采用的平均床面高程的變化量 |Δ(Z-h)|為一定時間間隔內(nèi)的沖淤厚度，總體上可反映出較長時段內(nèi)該斷面處于沖淤較小的狀態(tài)。經(jīng)統(tǒng)計，收集的枝城（2003—2011）、沙市（2003—2013，2015—2016）和監(jiān)利站（2003—2012，2015—2016）的數(shù)據(jù)共計3123 測次，符合判定標(biāo)準(zhǔn)的數(shù)據(jù)為561 測次。由于水位漲、落率較大的時段，河床變形一般會相對劇烈，故進(jìn)一步篩選了水位變幅的數(shù)據(jù)，共計453 測次。此外，由于選取的數(shù)據(jù)為枝城、沙市等水文站的水沙資料，其所在局部河段通常要求相對順直平整，水流集中，無整治工程，且河寬及水深等無明顯縱向變化（河流流量測驗規(guī)范（GB 50179-2015）），故可近似認(rèn)為本研究選取的資料均接近均勻流條件，一定程度上避免了水流非恒定性帶來的加速度對水位、水深等要素的影響，以及護(hù)灘護(hù)底等整治工程對床沙組成測量的影響。

然而，按照現(xiàn)行泥沙測量規(guī)范，一般采用五點法測量垂線平均含沙量，在距離河底0.5 m 近底范圍懸沙不采樣，這是目前含沙量測量廣泛存在的一個問題。He 等［17］等根據(jù)長江水文局荊江局在2010年后開展的荊江段近底含沙量測量數(shù)據(jù)，發(fā)現(xiàn)若忽略在近底區(qū)域（距離河底10%水深范圍內(nèi)）的含沙量，則泥沙通量在枝城、沙市和監(jiān)利站將分別被低估23.5%、9.4%和18.7%。因此，本研究計算的懸移床沙質(zhì)含沙量會較實際偏小?？紤]到在水沙數(shù)學(xué)模型應(yīng)用中，用于率定及驗證的含沙量同樣未考慮近底區(qū)泥沙。若要修正含沙量來率定挾沙力公式，則在模擬中的含沙量過程也需進(jìn)行相應(yīng)修正，否則仍無法準(zhǔn)確計算沖淤過程。鑒于上述情況，本文暫不進(jìn)行含沙量修正，有待于后續(xù)進(jìn)一步的研究。

3 張瑞瑾水流挾沙力公式參數(shù)確定

采用基于上述原則選取的相對沖淤平衡狀態(tài)下的含沙量資料（可近似等于水流挾沙力），對張瑞瑾挾沙力公式進(jìn)行率定，從而確定公式中的參數(shù)取值。主要計算步驟包括：懸移床沙質(zhì)含沙量S′與水沙綜合參數(shù)C′（扣除沖瀉質(zhì)部分）的計算；挾沙力公式中參數(shù)k和m與C′關(guān)系式的建立。

3.1 懸移床沙質(zhì)含沙量與水沙綜合參數(shù)的計算

（1）懸移床沙質(zhì)含沙量S′的計算。懸移床沙質(zhì)含沙量，即扣除沖瀉質(zhì)部分，為實際參與造床作用的含沙量S′=S-S沖瀉質(zhì)［15］。在具體計算中，通常將懸移質(zhì)級配曲線與相應(yīng)的床沙級配曲線進(jìn)行對比，來劃分懸移質(zhì)中的床沙質(zhì)與沖瀉質(zhì)部分［15］。具體劃分原則如下：在床沙級配曲線P＜10%的范圍內(nèi)，如出現(xiàn)比較明顯的拐點，就取與這一拐點相應(yīng)的床沙粒徑作為懸移質(zhì)泥沙中區(qū)分床沙質(zhì)與沖瀉質(zhì)的臨界粒徑dc。曲線中拐點的出現(xiàn)，表明懸移質(zhì)中大于此粒徑的泥沙是床沙中大量存在的，應(yīng)屬于床沙質(zhì)范圍；而小于此粒徑的泥沙是床沙中少有或沒有的，水流中這組泥沙幾乎不與床面發(fā)生交換，故屬于沖瀉質(zhì)范疇［15］。為簡化計算，通常取床沙級配曲線上P=5%對應(yīng)的粒徑作為臨界粒徑dc；然后dc在懸移質(zhì)級配曲線上對應(yīng)的百分比即為沖瀉質(zhì)所占比例，從而求出懸移床沙質(zhì)的含沙量。

采用該方法計算，第一組實測資料直接給出了三峽工程運用前長江中游10 個測站的懸移床沙質(zhì)含沙量S′，變化范圍在0.0026 ～0.4670 kg/m3之間（表1）。第二組數(shù)據(jù)中，懸移質(zhì)含沙量S則需根據(jù)上述原則扣除沖瀉質(zhì)部分：新廠站的床沙質(zhì)和沖瀉質(zhì)臨界粒徑dc范圍為0.010 ～0.104 mm；沖瀉質(zhì)所占比例ΔP沖為7% ～94%（平均值70%）；計算得到的懸移床沙質(zhì)含沙量S′范圍為0.0411 ～0.8334 kg/m3。第三組資料中，枝城、沙市、監(jiān)利3 站的懸移床沙質(zhì)含沙量S′范圍為0.0 ～0.3464 kg/m3；臨界粒徑dc的平均值分別為0.140、0.119、0.112 mm；且沖瀉質(zhì)所占比例分別介于25% ～99%（平均值82%）、18% ～96%（平均值57%）、12% ～93%（平均值47%）之間。

（2）水沙綜合參數(shù)C'的計算。計算水沙綜合參數(shù)時，需相應(yīng)扣除沖瀉質(zhì)部分，記為平均沉速ω′m采用韓其為方法計算，N為挾沙力分組數(shù)，ωk為第k粒徑組懸移質(zhì)泥沙沉速， ΔP?k為挾沙力級配，其值可由李義天［18］提出的方法確定。由于m是待確定量，故根據(jù)以往的研究成果將其分別取為0.5、1.0 及1.5，計算得到第三組的C′（453 測次）分別在0.024 ～13.775、0.024 ～12.253 及0.023 ～11.023 范圍內(nèi)?？梢姡S著m取值的增大，ω′m增大，C′ 相應(yīng)減小，但變幅不大，平均相對差值（MRE1.5-0.5）為9%。故此處將m取為1.0 進(jìn)行平均沉速的計算?？傮w上，長江中游的水沙綜合參數(shù)C′在10.0 以內(nèi)（表1），且其值在不同時段或不同位置差異較大。

3.2 挾沙力公式參數(shù)與水沙綜合參數(shù)的關(guān)系建立此處首先分析懸移質(zhì)總含沙量S與水沙綜合參數(shù)C（不扣除沖瀉質(zhì)）的關(guān)系，以此解釋選取水流挾沙力數(shù)據(jù)時需要扣除沖瀉質(zhì)部分的原因；然后進(jìn)一步研究懸移床沙質(zhì)含沙量S′與水沙綜合參數(shù)C′的關(guān)系；在此基礎(chǔ)上確定張瑞瑾水流挾沙力公式中參數(shù)k、m的計算關(guān)系。

（1）懸移質(zhì)含沙量S與水沙綜合參數(shù)C的關(guān)系。在河床處于相對沖淤平衡狀態(tài)且水流與河床中各組分泥沙充分交換時，認(rèn)為選取的懸移質(zhì)總含沙量S（不扣除沖瀉質(zhì)部分）即為水流挾沙力，那么挾沙力級配近似等于懸移質(zhì)級配，則根據(jù)懸移質(zhì)級配計算得到平均沉速ωm，進(jìn)而求得水沙綜合參數(shù)C。首先，點繪懸移質(zhì)總含沙量S與C的關(guān)系，如圖2（a）所示?？傮w上，數(shù)據(jù)點較為分散，擬合的冪函數(shù)關(guān)系的決定系數(shù)較低（R2=0.55）；此外，三峽工程運用后枝城和沙市站數(shù)據(jù)點與運用前兩站的數(shù)據(jù)點相比，更為散亂。主要是由于這里選取的河床沖淤幅度較小時的含沙量，是水流的實際挾沙力而不是理論挾沙力。水流的理論挾沙力為泥沙供應(yīng)充足情況下的水流挾沙力，該值僅由水流條件決定；實際挾沙力為在沖刷時河床補給受限條件下，水流實際能挾帶的含沙量，包括沖瀉質(zhì)和床沙質(zhì)部分；有效挾沙力則為真正參與造床作用的泥沙含量，為懸移質(zhì)中的床沙質(zhì)部分。此處在計算水沙綜合參數(shù)C時，假設(shè)了挾沙水流與河床中各組分泥沙是充分交換的，但實際上沖瀉質(zhì)泥沙無法從河床獲得補給，一直處于次飽和狀態(tài)，故選取沖淤平衡時的含沙量小于理論的水流挾沙力。

圖2 三峽工程運行前后長江中游水流挾沙力與水沙綜合參數(shù)的關(guān)系

由圖2（a）也可看出，枝城、沙市站的數(shù)據(jù)點較為散亂且偏低，主要是由于這兩站水流中沖瀉質(zhì)所占比例較大，實際挾沙力小于理論挾沙力。而三峽工程運用前各站和運行后監(jiān)利站的床沙組成較細(xì)，挾沙水流中沖瀉質(zhì)所占比例相對較小，故數(shù)據(jù)點較為集中?；谏鲜龇治?，選取含沙量和計算ωm時應(yīng)當(dāng)扣除沖瀉質(zhì)部分，即僅考慮真正與床沙進(jìn)行交換且使河床處于相對沖淤平衡狀態(tài)的泥沙（懸移質(zhì)中的床沙質(zhì)部分）。這樣便可排除計算ωm時考慮了沖瀉質(zhì)，但實際上沖瀉質(zhì)無法從河床得到補給而產(chǎn)生的誤差。上述即為計算水流挾沙力時要扣除沖瀉質(zhì)部分的原因。

（2）懸移床沙質(zhì)含沙量S′與水沙綜合參數(shù)C′的關(guān)系。通過上述分析可知，水流的有效挾沙力不僅取決于水流條件、懸沙級配，還受河床邊界條件（床沙組成）影響。河床邊界條件對水流挾沙力的影響，主要體現(xiàn)在對挾沙力級配ΔP?k的確定上。此處采用李義天方法［18］計算挾沙力級配，具體公式可寫為：

式中：ΔPbk為床沙級配；βk為與摩阻流速、 ΔPbk及ωk等相關(guān)的一個修正參數(shù)。采用式（2）計算挾沙力級配時，即認(rèn)為不考慮懸移質(zhì)中的沖瀉質(zhì)部分（ΔPbk在沖瀉質(zhì)范疇內(nèi)均為0），故沖瀉質(zhì)部分的水流挾沙力為0（S?k =S??ΔP?k）。而懸移質(zhì)含沙量S包含了沖瀉質(zhì)和床沙質(zhì)，故需將沖瀉質(zhì)扣除；并根據(jù)計算的挾沙力級配，計算平均沉速ω′m。

此處點繪了三峽工程運用前后長江中游懸移床沙質(zhì)含沙量S′與水沙綜合參數(shù)C′的關(guān)系（圖2（b）），可發(fā)現(xiàn)冪函數(shù)擬合曲線的決定系數(shù)（R2=0.76），較圖2（a）（R2=0.55）有顯著提高。直線斜率即為挾沙力公式的指數(shù)m≈1.5855，而通過擬合曲線與坐標(biāo)軸的交點則可求出系數(shù)k≈0.0454。這些結(jié)果較為合理，大致符合圖1的取值范圍。從時間尺度上看，三峽工程運用前，荊江段的有效水流挾沙力在0.0026 ～0.8334 kg/m3之間；而三峽工程運用后，S′則分布在0.0 ～0.3464 kg/m3之間，水流有效挾沙力總體上有所減小。從空間尺度上看，三峽工程運用后，有效水流挾沙力沿程遞增（S′枝城＜S′沙市＜S′監(jiān)利）。究其原因，主要是由于S′反映的是水流中有造床作用的床沙質(zhì)部分的挾沙力，三峽工程的運用使得河床發(fā)生粗化且越靠近大壩粗化程度越高，從而導(dǎo)致沖瀉質(zhì)臨界粒徑沿程增大，水流中床沙質(zhì)部分挾沙力相應(yīng)減小。

（3）挾沙力公式參數(shù)k和m計算關(guān)系的確定。從圖2（b）的點繪結(jié)果可看出，數(shù)據(jù)點的分布并非完全符合冪函數(shù)關(guān)系（直線），其變化趨勢總體呈現(xiàn)為斜率變小的曲線。故首先將S′ 與C′ 取對數(shù)，并進(jìn)行數(shù)據(jù)點平移，變形為lgS′+a和lgC′+b，旨在使兩個值均大于0；再在直角坐標(biāo)系中進(jìn)行擬合，并取最優(yōu)擬合方式；然后，求擬合曲線的斜率，即為不同水沙綜合參數(shù)C′ 下對應(yīng)的m值，并進(jìn)一步得到k值；最后，繪制k、m與C′ 的關(guān)系圖，并采用不同函數(shù)進(jìn)行分段擬合，取最優(yōu)擬合結(jié)果。根據(jù)上述方法，即可建立張瑞瑾挾沙力公式中參數(shù)k和m與C′ 之間的計算關(guān)系。此處a、b暫取為5 和2，首先采用各種函數(shù)形式進(jìn)行擬合，最終以對數(shù)擬合為優(yōu)（圖3），得到的關(guān)系式為1gS′+5=2.3097×ln(lgC′+2)+2.0641，并進(jìn)一步計算得到各水沙綜合參數(shù)C′ 對應(yīng)的k和m值。

圖3 三峽工程運用前后長江中游水流挾沙力（lgS ′+5）與水沙綜合參數(shù)（lgC ′+2）的關(guān)系

然后將數(shù)據(jù)點k、m與C′進(jìn)行分段擬合，使符合度達(dá)到最高（圖4），則不同C′范圍內(nèi)k、m的擬合曲線可分別表示為：

在長江中游河段，水沙綜合參數(shù)C′基本小于10.0，故獲得該范圍內(nèi)k、m的取值基本適用。當(dāng)C′＞10.0 時，其值可由張瑞瑾等［1］確定的關(guān)系曲線決定（圖1）。由圖4可知，系數(shù)k隨C′先減小后增大，當(dāng)0.2≤C′ ＜0.7 時，k由0.088 減小到0.044 kg/m3；而當(dāng)0.7≤C′ ≤10.0 時，k有逐漸增大趨勢，由0.046 增加至0.070 kg/m3。指數(shù)m隨C′的增大而減小，當(dāng)0.2≤C′≤10 時，m由1.86 減小至0.60?；谏鲜鍪剑?）和式（4），則可根據(jù)C′確定k、m值。

圖4 k、m 隨水沙綜合參數(shù)C ′的變化過程

此處采用長江中游各測站的低含沙量實測數(shù)據(jù)，確定了張瑞瑾水流挾沙力公式中參數(shù)k、m的計算關(guān)系，這些計算關(guān)系適用于低含沙水流及小的綜合水沙參數(shù)C′。該方法彌補了以往在低含沙條件下難以確定張瑞瑾挾沙力公式中參數(shù)的不足。但應(yīng)當(dāng)指出，這些計算關(guān)系僅適用于長江中游，在其他低含沙水流上不一定適用。由于水流的有效挾沙力與水沙及河床邊界條件均相關(guān)。故在不同的河流上，挾沙力公式的參數(shù)取值將有所不同。例如，韓其為挾沙力公式［19］中系數(shù)k取為0.245，偏大于本文k值，主要是由于其在公式率定時采用了長江、黃河、丹江口水庫等多條河流的資料，數(shù)據(jù)點范圍相對較廣，實際針對長江的k值應(yīng)顯著小于黃河上的值。然而，在邊界條件大致相同的情況下，這些計算關(guān)系仍具有參考價值。

4 張瑞瑾水流挾沙力公式計算實例對比

此處將張瑞瑾挾沙力公式中參數(shù)的計算式（3）和式（4），嵌入到已有的一維水沙數(shù)學(xué)模型中，進(jìn)一步驗證該公式的準(zhǔn)確性及適用性。該一維模型的水沙控制方程及數(shù)值解法，詳見文獻(xiàn)［20］。首先，按以往常規(guī)方法，將參數(shù)k和m取為常數(shù)，模擬了2016年荊江段（枝城至城陵磯）的水沙輸移過程；然后，由計算式確定參數(shù)k和m，再次模擬了該年荊江段水沙輸移過程；最后，比較了不同參數(shù)取值方法對模擬結(jié)果的影響。需說明，若研究河段受大規(guī)模護(hù)岸、護(hù)灘工程的守護(hù)而沖刷受限，可在模擬中通過修正挾沙力的方法來考慮這一影響［20］。

4.1 荊江河段概況荊江河段位于長江中游，上起枝城（荊3），下至城陵磯（荊186），全長347 km（圖5）。干流的水沙主要來自上游，其間通過松滋口、太平口、藕池口分流入洞庭湖，又于城陵磯處重新匯入長江干流。如忽略三口（松滋口、太平口及藕池口）分流的影響，枝城站可以作為控制斷面，則該站實測水沙過程代表進(jìn)入整個荊江段的水沙條件。三峽工程運用前（1994—2002年），枝城站年均水量為4304 億m3/a，而年均輸沙量為3.74 億t/a。工程運用后（2002—2017年），受人類活動及氣候變化的影響，枝城站年均水量減小到3952 億m3/a，為蓄水前的92%；受上游水土保持工程及三峽水庫攔沙作用的影響，枝城站年均輸沙量大幅減小到0.43億t/a，降幅高達(dá)88%［21］。此外，枝城站年均含沙量由工程運用前的0.869 kg/m3減小到運行后的0.109 kg/m3。三峽水庫的蓄水?dāng)r沙作用使得壩下游河床發(fā)生持續(xù)沖刷，2002—2017年荊江段平灘河槽累計沖刷量達(dá)10.51 億m3，單位河長的年均沖刷量為20.18 萬（m3/（a/km）），遠(yuǎn)大于蓄水前（1975—2002年）的3.18 萬（m3/（a/km））［21］。

圖5 荊江河段示意圖

在模型計算中，以枝城站日均的流量、含沙量和懸沙級配資料作為上游邊界條件；同時采用蓮花塘站的日均水位過程作為下游邊界條件；還需將三口分流的日均流量、含沙量和懸沙級配作為側(cè)向邊界條件。2016年枝城站最小、最大和平均流量分別為6950、34 000 和14 000 m3/s，且最小、最大和平均含沙量為0.003、0.147 及0.016 kg/m3；在蓮花塘站，水位在18.70～32.32 m 之間變化（圖6）。另外，采用荊江段173 個汛后實測固定斷面地形和72 個固定斷面的實測床沙級配（施測時間為2015年10月）作為初始河床邊界條件，其中床沙的中值粒徑范圍為0.17 ～0.33 mm；其余斷面的初始級配由這些實測值插值求得。由于在長江中游河段，推移質(zhì)輸沙量僅占懸移質(zhì)輸沙量的1%～2%［21-22］，對河床變形的影響很小，故在沖淤計算中未考慮推移質(zhì)輸移的影響。

4.2 k 和m 取常數(shù)時的計算結(jié)果分析荊江段2016年1月1日至12月31日模擬得到的各水文或水位站流量、水位過程與實測值十分吻合。在沙市站，實測最小、最大和平均流量分別為6030、28 600和12 610 m3/s，而相應(yīng)的計算值為6984、27 921 和12 802 m3/s；在監(jiān)利站，計算與實測的流量過程也幾乎一致，平均流量的相對誤差為3.9%。對于水位變化過程，在枝城、馬家店等水文或水位站計算與實測水位的平均相對誤差在0 ～1%之間。對于泥沙輸移計算，經(jīng)試算后確定當(dāng)張瑞瑾挾沙力公式參數(shù)k=0.05，m=1.55 時，計算得到的2016年荊江沙市及監(jiān)利站含沙量過程與實測過程符合相對較好，如圖7所示。沙市站最小、最大和平均含沙量實測值分別為0.011、0.328 和0.036 kg/m3，而相應(yīng)的計算含沙量為0.020、0.198 和0.039 kg/m3，且計算與實測含沙量的平均相對誤差（MRE）為49%；而在監(jiān)利站，計算與實測含沙量的平均相對誤差（MRE）為39%。此外，監(jiān)利站泥沙輸移情況較為特殊。2016年枝城、沙市及監(jiān)利站水位過程線的形狀基本一致；而對于含沙量，枝城站的含沙量過程線呈現(xiàn)枯水期小，汛期陡增的雙峰形狀（圖6），通常情況下，沙市及監(jiān)利站的含沙量過程應(yīng)與枝城站含沙量的變化趨勢一致，但2016年監(jiān)利站的含沙量在枯水期仍較大，且波動劇烈。這是由于監(jiān)利站的水力泥沙條件受洞庭湖出流頂托的影響顯著，含沙量變化較為復(fù)雜，故含沙量過程的模擬精度相對較低（圖7（b））。

圖7 2016年荊江段典型斷面含沙量計算與實測過程對比（k 和m 取常數(shù)）

4.3 k 和m 由公式確定時的計算結(jié)果分析此處將挾沙力公式的參數(shù)計算式，嵌入到了一維水沙數(shù)學(xué)模型中，使模型能根據(jù)各時刻計算得到的綜合水沙參數(shù)C' ，來確定參數(shù)k和m，從而求得各時刻的水流挾沙力。在此基礎(chǔ)上，重新模擬了2016年荊江段的水沙輸移過程?？芍?，沙市和監(jiān)利站計算與實測流量的平均相對誤差范圍為3% ～4%；而枝城、沙市及監(jiān)利站計算與實測水位的平均相對誤差均在0～1%之間（圖8）。因此，張瑞瑾挾沙力公式的參數(shù)取值采用不同方法時，對流量及水位計算結(jié)果的影響不大。然而沙市及監(jiān)利站計算與實測含沙量的平均相對誤差分別減小至26%和36%（圖9）。經(jīng)比較，計算含沙量與實測值的符合度較k和m取常數(shù)時基本一致或有所提高。因此，張瑞瑾挾沙力公式的參數(shù)取值采用公式計算時，能根據(jù)水力泥沙要素自動調(diào)整，對天然河流的適應(yīng)性更強，且含沙量過程的模擬也更為準(zhǔn)確。但總體上看，含沙量的計算和實測值誤差仍較大，尤其在監(jiān)利站，主要是因為該站的含沙量過程較為復(fù)雜，其較大程度上受到下游洞庭湖水沙入?yún)R的影響，僅根據(jù)現(xiàn)有的一維水沙數(shù)學(xué)模型無法考慮這一影響。

圖8 2016年荊江段典型斷面水位計算與實測過程對比

圖9 2016年荊江段含沙量過程計算與實測結(jié)果對比（k 和m 由公式?jīng)Q定）

分析2016年沙市、監(jiān)利站計算得到的綜合參數(shù)C′，其值范圍為0.26 ～0.89 和0.38 ～4.78。在這范圍內(nèi)，k的取值變幅不大，分別在0.044 ～0.070 kg/m3（平均值為0.052 kg/m3）和0.045～0.056 kg/m3（平均值為0.048 kg/m3）之間；而m變化較為顯著，范圍分別為1.21 ～1.72（平均值為1.46）和0.74 ～1.55（平均值為1.05）。通過比較，沙市站和監(jiān)利站的綜合水沙條件存在一定的差異，使得相應(yīng)的k和m取值也存在一定的差別。

5 結(jié)論

本研究基于長江中游的實測水沙資料，確定了低含沙量情況下張瑞瑾挾沙力公式中系數(shù)和指數(shù)的計算關(guān)系。并將該參數(shù)計算式嵌入到一維水沙動力學(xué)模型中，進(jìn)一步驗證公式的準(zhǔn)確性及適用性。本研究得到如下主要結(jié)論：

（1）選取三峽工程運行前、后相對沖淤平衡狀態(tài)下的水流含沙量（除去沖瀉質(zhì)部分）共計573組，將其近似等于水流挾沙力，并點繪水流挾沙力和綜合參數(shù)關(guān)系，從而確定參數(shù)k和m計算式。結(jié)果表明：系數(shù)k隨C′的增大呈現(xiàn)先減小后增大的趨勢，當(dāng)0.2≤C′＜0.7 時，其值由0.088 減小到0.044 kg/m3；而當(dāng)0.7≤C′ ≤10.0 時，該值由0.046 增加至0.070 kg/m3。指數(shù)m隨C′的增大而減小，在0.2≤C′≤10 范圍內(nèi)其值由1.86 減小至0.60。該參數(shù)計算式適用于低含沙水流及小的綜合水沙參數(shù)條件，彌補了以往在低含沙量條件下難以確定張瑞瑾挾沙力公式中參數(shù)的不足。

（2）將參數(shù)計算式應(yīng)用到一維水沙動力學(xué)模型中，并采用長江中游荊江段2016年的實測資料，對模型進(jìn)行驗證。結(jié)果表明：當(dāng)水流挾沙力公式的參數(shù)取為常數(shù)時，各水文站計算和實測含沙量平均相對誤差在39%～49%之間；而參數(shù)采用計算式確定時，含沙量平均相對誤差范圍為26%～36%。故挾沙力公式參數(shù)采用計算式確定時，相較于取為常數(shù)，模型能根據(jù)水力泥沙要素自動調(diào)整，對天然河流的適應(yīng)性更強，含沙量過程的模擬也更為準(zhǔn)確。