基于改進(jìn)馬斯京根法漓江上游河道流量演算模型探究

鄭金德，粟 忠，范玉潔，許景璇，李勇軍，蘇 釗，吳樹誠(chéng)，方榮杰

(1.桂林理工大學(xué)，廣西環(huán)境污染控制理論與技術(shù)重點(diǎn)實(shí)驗(yàn)室，廣西 桂林 541004；2.桂林理工大學(xué)，巖溶地區(qū)水污染控制與用水安全保障協(xié)同創(chuàng)新中心，廣西 桂林 541004；3.桂林水文中心，廣西 桂林 541001)

1 引言

桂林市以其獨(dú)特的人文和山水自然風(fēng)光而聞名于世，已被評(píng)為中外著名的國(guó)家級(jí)風(fēng)景和旅游城市。但是桂林市每年汛期的暴雨洪水幾乎都會(huì)給這座秀美的城市帶來(lái)不同程度的洪澇災(zāi)害，而且在枯水期漓江的水量又嚴(yán)重不足，給桂林市的城鄉(xiāng)飲用水、漓江旅游用水等造成一定的困擾。隨著旅游產(chǎn)業(yè)以及經(jīng)濟(jì)的發(fā)展，桂林市對(duì)防汛的要求也越來(lái)越高。桂林市急需準(zhǔn)確的雨情、水情信息，為桂林市的防汛工作提供相應(yīng)的參考依據(jù)。

目前，關(guān)于漓江流域河段流量演變規(guī)律的研究相對(duì)較少。馬斯京根法主要應(yīng)用于河段流量演算中，由McCarthy于1938年第一次提出。該方法從提出至今，一直被國(guó)內(nèi)外廣大從事水文工作的人員廣泛使用。馬斯京根模型的主要參數(shù)為K、x，對(duì)演算參數(shù)求解的方法大多采用人工試算法。在運(yùn)用馬斯京根模型進(jìn)行流量演算的過(guò)程中，發(fā)現(xiàn)演算出來(lái)的結(jié)果精度不是很高，為此許多科研工作者對(duì)馬斯京根法進(jìn)行了改進(jìn)。在對(duì)K、x求解的方面，JyS.Wu等提出了改進(jìn)圖解法， O，Donnell提出了一種直接三參數(shù)馬斯京根法，運(yùn)用矩陣計(jì)算方法直接推求K、x的值，且矩陣法較圖解法推求出來(lái)的結(jié)果更合理[1,2]。最小二乘法[3,4]是基于數(shù)學(xué)層面的計(jì)算方法，其概念清晰，在水文統(tǒng)計(jì)過(guò)程中被普遍應(yīng)用。例如，董曉華等[5]采用最小二乘法，對(duì)馬斯京根法河道流量演算精度的影響做了探究，并提出了采用最小二乘法有助于提高河道流量計(jì)算的精度；衛(wèi)小蓉[6]、薛夢(mèng)楠[7]等也基于最小二乘法對(duì)不同河道的流量進(jìn)行模擬，二者均取得了較好的演算結(jié)果。

本文針對(duì)漓江流域上游三街水文站至桂林水文站段干支流，選取三街水文站、靈川水文站和桂林水文站的逐日實(shí)測(cè)流量數(shù)據(jù)，在有約束條件下推導(dǎo)出流量演進(jìn)系數(shù)的最優(yōu)估計(jì)值，客觀地建立適用于漓江流域上游的馬斯京根流量演算模型，并探究該河段的流量演進(jìn)規(guī)律，進(jìn)而為漓江流域防洪工作地進(jìn)行及水資源的合理調(diào)度提供相應(yīng)的理論依據(jù)。

2 材料與方法

2.1 流域概況

漓江流域上游桂林水文站以上的集水區(qū)總面積為2762 km2，平均海拔約150 m，河段主要干流長(zhǎng)105 km，流經(jīng)桂林市區(qū)的河段長(zhǎng)4.93 km，整個(gè)流域的地勢(shì)自西北向東南傾斜[8]。且桂林水文站以上流域內(nèi)雨水充足，水資源豐富；流域多呈羽狀，支流短，河道類型為峽谷式河槽，河床的坡降大，漓江流域上游存在多個(gè)暴雨區(qū)。桂林市氣候?qū)賮啛釒Ъ撅L(fēng)氣候，多年平均降雨量為1889.4 mm，汛期大多集中在3～8 月份[9]，汛期降雨量約占年降水量的75%，在汛期內(nèi)漓江流域降水較多，頻繁發(fā)生水澇災(zāi)害。桂林市區(qū)三街水文站至桂林水文站主要河流基本特征見表1，河流及水文站網(wǎng)分布見圖1。

表1 漓江桂林水文站以上主要河流基本特征

2.2 數(shù)據(jù)來(lái)源

本次模型探究所使用的資料由桂林水文中心所提供，主要研究河段為漓江上游的三街水文站至桂林水文站，該河段內(nèi)甘棠江為漓江的一級(jí)支流，甘棠江上布設(shè)有靈川水文站，其它支流并未布設(shè)有相應(yīng)的水文測(cè)站。三街水文站處于漓江上游河段的干流上，該水文站以上的流域內(nèi)建有斧子口水庫(kù)、川江水庫(kù)和小溶江水庫(kù)，各水庫(kù)所處位置見圖1，且流域以上支流眾多，地形地勢(shì)復(fù)雜。進(jìn)行流量演進(jìn)的前提是演算河道上下游的流量資料要相互對(duì)應(yīng)，三街水文站為新建的水文站，所具有的流量監(jiān)測(cè)數(shù)據(jù)年份較少，并且由前文可知桂林市的汛期主要集中在每年的3～8月份。因此，本次研究所使用的數(shù)據(jù)為三街水文站、靈川水文站和桂林水文站2019、2020年3～8月份的逐日實(shí)測(cè)流量數(shù)據(jù)。以2019年3～8月份的逐日實(shí)測(cè)流量數(shù)據(jù)進(jìn)行模型參數(shù)的率定，以2020年3～8月份的逐日實(shí)測(cè)流量數(shù)據(jù)對(duì)模型進(jìn)行了驗(yàn)證。

圖1 桂林市區(qū)河流及水文站網(wǎng)分布

2.3 馬斯京根流量演算模型

馬斯京根模型主要是一種在河道流量(洪水)演算中被廣泛應(yīng)用的模型，它利用水量平衡方程式代替連續(xù)方程，以槽蓄方程代替復(fù)雜的水動(dòng)力方程，極大地簡(jiǎn)化了計(jì)算過(guò)程[10]。本文依據(jù)漓江流域上游河段干支流洪水匯集的特性，采用先合后演[11]的流量演算方法，將干、支流斷面上三街水文站與靈川水文站的實(shí)測(cè)流量相互疊加在一起，以此為基礎(chǔ)作為漓江流域上游河段的總?cè)肓髁浚凑諢o(wú)支流河段的演算方法進(jìn)行連續(xù)馬斯京根演算。

2.3.1 馬斯京根方程

在第1、2時(shí)段內(nèi)對(duì)水量平衡方程和槽蓄方程[12,13]進(jìn)行差分求解，可推導(dǎo)出流量演算方程式：

Q2=C0I2+C1I1+C2Q1

(1)

式(1)中：

(2)

(3)

(4)

C0+C1+C2=1.0

(5)

因此，在使用馬斯京根方程時(shí)，只要確定了參數(shù)K、x與△t的值，就可根據(jù)式(2)、(3)、(4)分別推導(dǎo)得出C0、C1、C2的值，然后再根據(jù)流量演算方程式(1)演算出下斷面的最終流量過(guò)程Q(t)。

顯然，確定模型演算參數(shù)x、K的值，成為使用馬斯京根流量演算模型的關(guān)鍵。模型參數(shù)x、K值的確定，一般根據(jù)馬斯京根法W～Q`呈單一直線的基本假設(shè)對(duì)x、K的值進(jìn)行確定，傳統(tǒng)的確定x、K值的方法有試錯(cuò)法、最小二乘法、最小面積法等。利用傳統(tǒng)的方法推求x、K值是使W～Q`呈單一直線的最優(yōu)估計(jì)值，但是根據(jù)所推導(dǎo)出的x、K值計(jì)算的流量演算系數(shù)未必能使演算出流量與實(shí)測(cè)出流量的擬合誤差最小[14]。而基于最小二乘法，直接對(duì)流量演算系數(shù)C0、C1、C2的最優(yōu)估值進(jìn)行推導(dǎo)，最后再根據(jù)相應(yīng)公式進(jìn)行x、K值的反算，再利用所推求出的系數(shù)進(jìn)行河道流量演算，則可取得較好的演算結(jié)果。

2.3.2 模型參數(shù)優(yōu)選依據(jù)

模型參數(shù)的優(yōu)選以演算出流量與實(shí)測(cè)出流量的離差平方和最小為依據(jù)，構(gòu)造目標(biāo)函數(shù)[15]表達(dá)式(6)與約束方程表達(dá)式(7)：

(6)

(7)

式(6)中S為離差平方和，n為實(shí)測(cè)流量點(diǎn)數(shù)。

由式(2)、(3)、(4)可知C0、C1、C2必須滿足式(5)的要求。因此可以根據(jù)拉格朗日乘數(shù)法構(gòu)造一個(gè)非線性規(guī)劃問(wèn)題，由式(6)、式(7)構(gòu)成等式約束方程，再構(gòu)造拉格朗日函數(shù)，將等式所需要約束的極值問(wèn)題轉(zhuǎn)化為一個(gè)無(wú)線性約束的極值問(wèn)題，其極值函數(shù)為式(8)：

L(λ,C0,C1,C2)=

(8)

式(8)中，λ為拉格朗日乘數(shù)。式(8)達(dá)到最小的C0、C1、C2值即為式(7)的解，根據(jù)極值存在的必要條件，則有式(9)、(10)：

(9)

(10)

經(jīng)處理化簡(jiǎn)后得式(11)、(12)、(13)：

(11)

(12)

(13)

聯(lián)立式(5)與式(11)、(12)、(13)即可求出流量演算系數(shù)C0、C1、C2的最優(yōu)估值。然后根據(jù)式(14)、(15)即可反求出x、K的值。

(14)

(15)

2.4 模型的驗(yàn)證方法

本文主要以確定性系數(shù)[16]和相關(guān)系數(shù)[17]作為模型預(yù)測(cè)精度的評(píng)價(jià)標(biāo)準(zhǔn)，確定性系數(shù)DC按式(16)計(jì)算，相關(guān)系數(shù)R按式(17)計(jì)算：

(16)

(17)

3 結(jié)果與分析

3.1 馬斯京根流量演算系數(shù)的推求

選取三街水文站、靈川水文站至桂林水文站2019年3～8月份的逐日實(shí)測(cè)流量資料，代入式(8)對(duì)馬斯京根模型的流量演算系數(shù)C0、C1、C2進(jìn)行推求，然后根據(jù)式(14)與式(15)推求出K、x的值；在對(duì)傳統(tǒng)的馬斯京根模型的流量演算系數(shù)進(jìn)行推求，兩者所推求出的流量演算系數(shù)值見表2。把推求得到的C0、C1、C2值代入式(1)進(jìn)行模型率定期的流量演算，演算出來(lái)的流量值與實(shí)測(cè)流量值成果見圖2。

表2 河段流量演算系數(shù)值

圖2 演算流量與實(shí)測(cè)流量對(duì)比(率定期)

由圖2可知，利用上文所介紹的方法與傳統(tǒng)的方法對(duì)馬斯京根法的流量演算系數(shù)C0、C1、C2進(jìn)行最優(yōu)取值，然后在進(jìn)行河道流量演算，得出實(shí)測(cè)流量值與演算流量值具有大致相同的變化趨勢(shì)，演算出的流量值與水文站實(shí)測(cè)的流量值擬合較好。

3.2 模型檢驗(yàn)

把由2019年3～8月份的逐日實(shí)測(cè)流量數(shù)據(jù)推求出的河段流量演算系數(shù)值，代入到三街水文站至桂林水文2020年3～8月份的逐日流量數(shù)據(jù)進(jìn)行模擬，模擬結(jié)果見圖3。從圖2、圖3可以看出，無(wú)論是在模型的率定期或者驗(yàn)證期，所演算得出的流量值與實(shí)際監(jiān)測(cè)得到的流量值之間是具有較大偏差的。在流量較小的階段，演算值與實(shí)測(cè)值之間可以達(dá)到很好的擬合效果，而在流量較大的階段，演算值與實(shí)測(cè)值之間盡管存在著大體一致的變化趨勢(shì)，但是仍然存在著一定的偏差。產(chǎn)生誤差的主要原因是，在進(jìn)行流量演算的過(guò)程中并未加入桃花江與南溪河的流量進(jìn)行運(yùn)算，而由表1 可知桃花江的多年平均流量為12.2 m3/s，如若加上桃花江與南溪河的流量進(jìn)行演算，則可取得更好的結(jié)果。

由圖3和表3可以看出，把由2019 年的流量數(shù)據(jù)確定出來(lái)的馬斯京根演算系數(shù)C0、C1、C2值，代入到三街水文站2020 年的流量監(jiān)測(cè)數(shù)據(jù)進(jìn)行演算，演算出的結(jié)果與桂林水文站實(shí)測(cè)的流量結(jié)果擬合程度較好，二者都達(dá)到了乙級(jí)標(biāo)準(zhǔn)，可為洪水預(yù)報(bào)工作提供相應(yīng)的參考依據(jù)。

由表3可以明顯看出，基于改進(jìn)的馬斯京根法所推求出來(lái)的流量值，不論是在率定期還是驗(yàn)證期，其確定性系數(shù)與相關(guān)系數(shù)都比由馬斯京根法推求出來(lái)的高，但是對(duì)比二者的平均相對(duì)誤差值，改進(jìn)的馬斯京根法卻比馬斯京根法的高。從圖2、圖3可以看出不論是在模型的率定期還是驗(yàn)證期，改進(jìn)的馬斯京根法與實(shí)測(cè)流量值之間具有更好的擬合效果，曲線的變化趨勢(shì)更加一致。改進(jìn)的馬斯京根法是基于最小二乘法，直接推導(dǎo)出馬斯京根模型的河道流量演算系數(shù)，其計(jì)算過(guò)程簡(jiǎn)單，取得的擬合效果更好。而由馬斯京根法推求的流量演算系數(shù)，計(jì)算過(guò)程繁瑣，演算系數(shù)的推求受計(jì)算者的影響較大，所取得的演算結(jié)果往往與實(shí)測(cè)的結(jié)果擬合效果較差。因此，改進(jìn)的馬斯京根法比馬斯京根法更適合應(yīng)用于漓江上游的流量演算。

圖3 演算流量與實(shí)測(cè)流量對(duì)比(驗(yàn)證期)

表3 模型精度評(píng)定結(jié)果

4 結(jié)論

(1)基于改進(jìn)的馬斯京根流量演算法，推求出了適合漓江上游流域的流量演算系數(shù)C0、C1、C2(0.37；0.62；0.01)的值。改進(jìn)的馬斯京根流量演算法較傳統(tǒng)的馬斯京根法取得更好的擬合效果，在模型驗(yàn)證期內(nèi)，改進(jìn)的馬斯京根法所推求出的結(jié)果確定性系數(shù)達(dá)到了0.812，相關(guān)系數(shù)達(dá)到了0.940，平均相對(duì)誤差為-0.195，達(dá)到了乙級(jí)精度。可為桂林市的洪水預(yù)報(bào)作業(yè)提供相應(yīng)的參考依據(jù)。

(2)改進(jìn)的馬斯京根法較傳統(tǒng)的馬斯京根法，減少了繁瑣的參數(shù)試算過(guò)程，提高了參數(shù)的選取效率，同時(shí)也增強(qiáng)了演算結(jié)果的可靠性。

(3)本次模型探究所使用的數(shù)據(jù)是水文測(cè)站每日監(jiān)測(cè)的流量數(shù)據(jù)，其精度可能會(huì)存在一定的問(wèn)題，且使用的數(shù)據(jù)資料年份較少，會(huì)給參數(shù)的確定帶來(lái)一定的影響，再加上有一些支流上并未布設(shè)水文測(cè)站，所以在今后的工作中應(yīng)及時(shí)掌握流域的來(lái)水情況，及時(shí)地對(duì)模型參數(shù)進(jìn)行相應(yīng)的修正。