王船海，俞 悅，吳金寧，陳 鋼，馬騰飛，曾賢敏

(1.河海大學 水文水資源學院，江蘇 南京 210098；2.河海大學 水文水資源國家重點實驗室，江蘇 南京 210098；3.江蘇省水文水資源勘測局常州分局，江蘇 常州 213022)

1 研究背景

河道匯流演算屬于水力學問題[1]，但由于對河道實測斷面資料的要求高、計算耗時長等問題，水文學家們研究和發(fā)展了一系列的水文學方法和簡化的水動力學方法來進行演算[2-3]。馬斯京根法便是最經(jīng)典的水文學方法之一。水動力學方法[4-6]則主要是通過求解圣維南方程組，但由于圣維南方程組屬于復雜的偏微分方程，求解過程較為復雜[7-8]。從首次面世以來，馬斯京根法經(jīng)歷了從經(jīng)驗性方法上升為具有物理基礎方法的發(fā)展歷程[9]。在長期發(fā)展的過程中，研究者們在傳統(tǒng)線性方法的基礎上，發(fā)展出一系列的非線性方法，如馬斯京根-康吉法[10]、可變參數(shù)馬斯京根-康吉法[11]等。馬斯京根法計算過程簡便，所需資料少，且精度能滿足一般科研、工程需求。目前，馬斯京根法在我國已經(jīng)取得了廣泛應用，利用馬斯京根法進行匯流演算的成果非常豐富[12-14]。同時，經(jīng)過廣大學者的不斷探索，對圣維南方程組進行求解的水動力學模型發(fā)展已經(jīng)日趨成熟，能夠快速有效地模擬河道流量、水位、流速等要素。

“等效河道”是一種將復雜的實際河流簡化為具有相似水動力特性的“理想河道”的概念，能夠將復雜的自然河道形態(tài)簡化為等效的理想河道形狀和參數(shù)，與原河道具有相似的水動力特性的虛擬河道。早在1988年Chow等[15]就已提出了等效河道的概念，然而由于河道水力特性復雜，等效的河道需要考慮的要素過多。目前國內(nèi)外對等效河道的研究較少，高學平等[16]通過在河道內(nèi)分區(qū)設置植物，對含植物河道等效創(chuàng)面阻力進行了試驗性研究；龔定等[1]根據(jù)原河道大斷面形狀對河道斷面進行了梯形概化?？梢钥闯?，關于等效河道的研究大部分是采取的試驗性研究，或是對河道斷面直接進行簡單概化處理，而通過簡便的參數(shù)實現(xiàn)等效河道的公式推導的研究基本沒有。

隨著國家數(shù)字孿生流域建設中“四預”的新要求，對洪水預報提出了新的挑戰(zhàn)，包括預報河道水位、流速等信息。目前，對馬斯京根法的相關研究大多聚集在對K和X參數(shù)值更精確的獲取[17-18]，以及馬斯京根法演算結果精度的進一步提高中[19-21]。傳統(tǒng)的馬斯京根法輸出結果類型單一，對河道內(nèi)其他要素的獲取需更多資料支撐，無法與水動力模型相結合進行運用。因此，本文提出了一種基于馬斯京根參數(shù)K、X的等效河道計算方法，通過構建等效河道，可求解出資料缺乏地區(qū)河道的河寬、河道比降、最大水深信息，通過水動力模型進行演算后實現(xiàn)了河道水位、流速、流量等要素的模擬。

2 計算原理

2.1 傳統(tǒng)馬斯京根法水量平衡方程與槽蓄方程聯(lián)立：

(1)

得：

Q2=C0·I2+C1·I1+C2·Q1

(2)

其中：

式中：I0為上斷面流量，I1、I2分別為時段始末上斷面流量，m3/s；Q0為下斷面流量，Q1、Q2分別為時段始末下斷面流量，m3/s；W0為河段槽蓄量，W1、W2分別為時段始末河段槽蓄量，m3；Δt為計算時段；K和X為馬斯京根法主要參數(shù)。K為蓄量常數(shù)，表示穩(wěn)定流時河段的傳播時間，隨流量的變化而變化；X為楔蓄系數(shù)，由河道楔蓄因素和調(diào)蓄因素組成。

X的計算公式如下：

(3)

式中：l為特征河長，m；L為總河段長，m。

2.2 總河段K總、X總與分河段Ki、Xi關系如圖1所示，將一個長河段分為N個子河段，建立馬斯京根法總河段K總、X總與分河段Ki、Xi關系。

圖1 馬斯京根分河段示意

(4)

(5)

假設每個分河段的K1=K2=…=KN，X1=X2=…=XN，則：

(6)

式中：K總、X總為馬斯京根總河段參數(shù)；Ki、Xi為馬斯京根分河段參數(shù)。

2.3Ki、Xi與水動力模型河道要素轉換關系結合特征河長，利用馬斯京根法分河段參數(shù)Ki、Xi建立等效河道要素的轉換關系。

已知曼寧公式：

(7)

波速-河長傳播時間關系式：

Li=Cs·Ki=Ki·αu

(8)

特征河長計算公式：

(9)

令水力半徑R=R(H)，聯(lián)立式(7)(8)，得到水力半徑R、河道比降i0與馬斯京根法分河段參數(shù)Ki之間的關系：

(10)

斷面面積A=f(H)。將曼寧公式與流量公式(Q=f(H)·u)聯(lián)立，得洪峰流量Q的表達式：

(11)

對式(11)中的H偏導得：

(12)

通過聯(lián)立式(3)和式(9)，構建等效河道。

將式(12)代入式(9)，得：

(13)

聯(lián)立式(3)與式(13)，得：

(14)

建立河道比降i0與馬斯京根參數(shù)之間的關系。由式(11)得到以下變形：

(15)

將式(15)代入式(14)得：

(16)

由式(10)、式(11)、式(16)得到通用斷面形狀的馬斯京根法參數(shù)Ki、Xi與等效河道要素之間的轉換關系。

3 具體河道斷面計算方法

3.1 矩形斷面如圖2所示，假設河道斷面形狀為矩形時，經(jīng)簡化，水力半徑Rr=Hr，斷面面積f(H)r=Br·Hr，f′(H)r=Br，代入式(16)得河道比降i0r：

i0r=Dr·Hr

(17)

圖2 矩形斷面示意

式中：

(18)

將i0r=Dr·Hr帶入式(10)，得斷面最大水深Hr：

(19)

由式(11)得河寬Br：

(20)

由式(17)(19)(20)得到矩形斷面的馬斯京根法分河段參數(shù)Ki、Xi與水動力模型河寬Br、河道比降i0r、最大水深Hr之間的轉換關系。

圖3 拋物線型斷面示意

(21)

代入式(16)得河道比降i0p：

i0p=Dp·Hp

(22)

式中：

(23)

將i0p=Dp·Hp帶入式(10)，得斷面最大水深Hp：

(24)

由式(11)得β的計算公式：

(25)

則河寬Bp：

(26)

由式(22)(24)(26)得到拋物線型斷面的馬斯京根法分河段參數(shù)Ki、Xi與水動力模型河寬Bp、河道比降i0p、最大水深Hp之間的轉換關系。

根據(jù)Ki、Xi與等效河道要素的轉換關系，通過馬斯京根參數(shù)Ki、Xi、總河段長L、洪峰流量Q、糙率n、平均流速轉換波速系數(shù)α及分河段數(shù)N，可以計算得到不同概化斷面下等效河道的河寬B、河道比降i0、最大水深H等河道信息。本文只提出了拋物線型和矩形這兩種概化斷面下的轉換關系，根據(jù)相同的原理，其他形狀的概化斷面(如三角形、梯形等)也可推導求出。

4 基于K、X等效河道的具體應用及合理性分析

4.1 基于K、X等效河道的具體應用何惠等[22]提出了用最小二乘法對馬斯京根法的K、X參數(shù)進行最優(yōu)估計，并以1961年海河流域南運河稱溝灣至臨清段的一次洪水過程為例，其中K=13.05、X=-0.2716?，F(xiàn)以該方法率定的K、X值運用于本文提出的等效河道方法，對該場次洪水進行演算。

南運河位于河北省南部，南運河稱溝灣至臨清段全長83.8 km，如圖4(a)所示。本文提出的等效河道方法對該河段的概化如圖4(b)(c)所示。整個河道概化為一條長83.8 km的長河段，取N為100，將河段等距離劃分為100個子河段，取糙率n為0.02，洪峰流量Q為550 m3/s。經(jīng)計算，當河道斷面概化為矩形時，河道全程比降為0.0045%，河寬為36.04 m；當河道斷面概化為拋物線型時，河道全程比降為0.0042%，河寬為51.62 m。K、X值對應的等效河道信息，轉換后可用于水動力模型的邊界條件。由臨清斷面形狀信息，結合曼寧公式及流速流量關系可推出臨清斷面水位流量關系。取南運河稱溝灣實測流量為上邊界條件，臨清斷面的水位流量關系為下邊界條件進行模擬。

圖4 實際河道及下邊界水位流量關系

圖5和表1為出口臨清斷面實測流量、參考文獻[22]方法演算流量及本文方法演算流量對比。從表1統(tǒng)計結果來看，矩形斷面下，均方根誤差RMSE為10.71，比文獻方法大4.64；可決系數(shù)R2為0.99，與文獻方法相近，均在0.99 以上。拋物線型斷面下，均方根誤差RMSE為10.76，比文獻方法大4.69；可決系數(shù)R2為0.99，與文獻方法相近，均在0.99 以上。本算例下，兩種概化方式演算結果相差不大，矩形的概化方式表現(xiàn)更好。從圖5流量過程線來看，本文提出的兩種概化方法演算的臨清流量結果較好，與文獻[22]方法演算結果相差不大，能反應臨清出流過程。

圖5 稱溝灣—臨清流量過程線(1961年)

表1 秤鉤灣—臨清流量計算結果對比

本文方法演算出不同斷面處的流速情況對比如圖6，可以看到，兩種概化方式下，臨清斷面處流速情況與圖5中的流量過程線一致，計算結果具有可靠性。由圖6中不同斷面的流速情況可以看出，矩形斷面下，各斷面流速于8月21日達到最大，其中臨清斷面流速為1.6 m/s。同時，拋物線型斷面下，各斷面在8月22日達到最大流速，其中臨清處流速最大，為1.37 m/s。

圖6 不同等效河道斷面處流速情況

4.2 合理性分析本文提出的方法使用的是洪峰流量，而非入流過程，不同流量下相應的河道概化斷面也會出現(xiàn)偏差，選用洪峰流量進行演算，是考慮到洪峰流量是最值得關注的特殊情況。因此，不同入流及時間變化都是非線性分析的重要考慮因素。從表2和圖7可知，當輸入的洪峰流量分別為200，400，600 m3/s，概化斷面為拋物線型時，3種擬合效果均較好，R2均趨于1，且RMSE均小于16。但不同的洪峰流量在計算得到的洪水漲落過程也呈現(xiàn)出一定的差異性，體現(xiàn)了非線性的特征。從這三種計算結果來看，當洪峰流量為400 m3/s時，誤差最小。

圖7 不同洪峰流量對應的模擬結果

表2 不同洪峰流量Q下計算流量結果對比

為了進一步驗證本文方法計算水位、流速等信息的準確及合理性，將本文提出的兩種概化河道計算出的結果與臨清斷面實測日平均水位進行對比。如表3所示，矩形斷面下，均方根誤差RMSE為0.90，可決系數(shù)R2為0.97；拋物線型斷面下，均方根誤差RMSE為0.95，可決系數(shù)R2為0.95。如圖8所示，本文提出的兩種概化方式均能在一定程度上反應水位的變化情況，趨勢一致。

圖8 臨清斷面實測水位與計算水位對比

表3 臨清斷面水位計算結果對比

計算不同洪峰流量下，兩種河道斷面概化方式對水位的計算結果如表4和圖9所示?？梢钥吹?，不同的洪峰流量對水位的計算結果呈現(xiàn)出一定的非線性特征，但從RMSE和R2這兩個指標來看，模擬結果表現(xiàn)優(yōu)秀，RMSE均小于1.4，R2均趨于1。從水位過程線中可以清晰的看出，不同洪峰流量下，演算得出的結果差異較小?？偟膩砜?，相較于拋物線型斷面，矩形斷面計算結果更優(yōu)。

圖9 不同洪峰流量下臨清斷面實測水位與計算水位過程

表4 不同洪峰流量下臨清斷面實測水位與計算水位結果對比

實測的河道大斷面資料如圖10所示，可以看到，稱溝灣測站過水斷面形狀較不規(guī)則，主槽間距為52.78 m。在最大水位為41.18 m時，本文方法概化的矩形河道斷面河寬36.04 m，拋物線型河道斷面河寬44.67 m?？梢钥吹剑@兩種概化方式都十分對稱，是一種理想化的河道概化方式，且大致能夠描述河道斷面。值得一提的是，本文提出的概化河道與實測大斷面仍存在一定差距，這主要由于本文提出的概化方法使整個河道每個斷面都均化為唯一的斷面形狀，稱溝灣斷面只是其中一個斷面，因此存在一定差距是合理的現(xiàn)象。

圖10 稱溝灣站實測大斷面形狀與概化斷面對比

5 結論

本文提出了基于K、X等效河道的計算方法，通過馬斯京根法中參數(shù)K、X構建等效河道，推演出不同的概化河道斷面(矩形、拋物線型)下的河寬、河道比降和最大水深計算公式，并根據(jù)參考文獻中率定好的參數(shù)在實例中進行驗證。結果表明，該方法在傳統(tǒng)馬斯京根法的基礎上，通過上斷面流量、馬法參數(shù)K、X值、河道長度L、糙率n、平均流速轉換波速系數(shù)α和洪峰流量Q，就可以推演出概化河道斷面的河寬B、河道比降i0和最大水深H推導公式。利用這些信息，除了可以通過水動力學的方法對資料缺乏地區(qū)的河道進行流量預報外，在進一步的應用中，任意時刻、河道任意斷面水深、流速信息也能在水動力模型中模擬出來。總體來說，本研究實現(xiàn)了河道大斷面資料的擴充，改進了傳統(tǒng)的水文學河道匯流計算模式，有利于實現(xiàn)河道水位、流速過程模擬，為傳統(tǒng)的馬斯京根法轉換為水動力模型計算提供了新方法，具有物理意義明確、計算方法簡便、計算精度高等優(yōu)點。

值得注意的是，傳統(tǒng)的馬斯京根法的分河段是為了解決長河段下出現(xiàn)的非線性的問題，而本文提出的分河段則是轉換后用于水動力模型計算的，N可以為任意值。本研究主要提出的是線性馬斯京根法下的等效河道概化方法，但由于實際河道水力要素復雜，往往都會或多或少的受到非線性因素的影響，K、X參數(shù)的值會隨著Q的變化而變化，因此非線性馬斯京根法轉換為等效河道的方法也值得研究，以進一步減小計算誤差，并縮小概化斷面形狀與實測斷面形狀之間的差距。