基于馬斯京根參數(shù)K、X的等效河道計(jì)算方法研究

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

(1.河海大學(xué) 水文水資源學(xué)院，江蘇 南京 210098；2.河海大學(xué) 水文水資源國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210098；3.江蘇省水文水資源勘測(cè)局常州分局，江蘇 常州 213022)

1 研究背景

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

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

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

2 計(jì)算原理

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

(1)

得：

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

(2)

其中：

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

X的計(jì)算公式如下：

(3)

式中：l為特征河長(zhǎng)，m；L為總河段長(zhǎng)，m。

2.2 總河段K總、X總與分河段Ki、Xi關(guān)系如圖1所示，將一個(gè)長(zhǎng)河段分為N個(gè)子河段，建立馬斯京根法總河段K總、X總與分河段Ki、Xi關(guān)系。

圖1 馬斯京根分河段示意

(4)

(5)

假設(shè)每個(gè)分河段的K1=K2=…=KN，X1=X2=…=XN，則：

(6)

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

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

已知曼寧公式：

(7)

波速-河長(zhǎng)傳播時(shí)間關(guān)系式：

Li=Cs·Ki=Ki·αu

(8)

特征河長(zhǎng)計(jì)算公式：

(9)

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

(10)

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

(11)

對(duì)式(11)中的H偏導(dǎo)得：

(12)

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

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

(13)

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

(14)

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

(15)

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

(16)

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

3 具體河道斷面計(jì)算方法

3.1 矩形斷面如圖2所示，假設(shè)河道斷面形狀為矩形時(shí)，經(jīng)簡(jiǎn)化，水力半徑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與水動(dòng)力模型河寬Br、河道比降i0r、最大水深Hr之間的轉(zhuǎn)換關(guān)系。

圖3 拋物線型斷面示意

(21)

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

i0p=Dp·Hp

(22)

式中：

(23)

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

(24)

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

(25)

則河寬Bp：

(26)

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

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

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

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

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

圖4 實(shí)際河道及下邊界水位流量關(guān)系

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

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

表1 秤鉤灣—臨清流量計(jì)算結(jié)果對(duì)比

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

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

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

圖7 不同洪峰流量對(duì)應(yīng)的模擬結(jié)果

表2 不同洪峰流量Q下計(jì)算流量結(jié)果對(duì)比

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

圖8 臨清斷面實(shí)測(cè)水位與計(jì)算水位對(duì)比

表3 臨清斷面水位計(jì)算結(jié)果對(duì)比

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

圖9 不同洪峰流量下臨清斷面實(shí)測(cè)水位與計(jì)算水位過(guò)程

表4 不同洪峰流量下臨清斷面實(shí)測(cè)水位與計(jì)算水位結(jié)果對(duì)比

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

圖10 稱溝灣站實(shí)測(cè)大斷面形狀與概化斷面對(duì)比

5 結(jié)論

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

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